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Abstract 


SERDP project MM-1537 entitled “Handheld Frequency Domain Vector EMI Sensing for UXO Dis- 
crimination” is complete. The original objectives of MM1537 consisted of both hardware and software 
goals. On the hardware side, the goal was to develop an innovative vector (multi-axis) handheld UWB 
electromagnetic induction sensor in the frequency domain, with precise 3-D positioning, for close inter- 
rogation of anomalies. This new instrument would allow a new higher level of UXO discrimination in 
the vicinity of a noteworthy magnetic response. Data processing software accompanying the new instru- 
ment would also be developed. The goal here was to develop clutter-tolerant signal processing for UXO 
discrimination using the data provided by the new sensor, based on new, high-fidelity, physically com- 
plete forward modeling (the Standardized Excitation Approach), rigorous instrument characterization, 
and on new processing techniques. 

We have constructed the GEM-39* along with a positioning system with subcentimeter accuracy at 
a range of about 2 meters. This new instrument has been successfully used to acquire data both under lab 
conditions and in test plots with emplaced targets. Results from all blind tests (see Sec. 5) indicate the 
GEM-39* can acquire diverse and accurate vector data with a SNR similar to or better than prior GEM 
models. Single target inversion results for laboratory blind tests were 100% accurate. Blind tests from 
the test plots at CRREL, with either one or two targets in close proximity, were also very encouraging 
though less accurate. 

A novel “beacon” positioning, which uses the primary field of the GEM-39* itself to locate the 
sensor head, was developed and successfully deployed. Early ground and stiffness sensitivities were 
overcome with limited redesigning of the main sensor head. The GEM-39* successfully acquires three 
axis data of the secondary field near the center of the instrument head while remaining portable for 
operators and reproducable by Geophex Ltd.. Details of the GEM-3Y* hardware development can be 
found in Sec. 3. 

Our group has developed algorithms and models which accompanied the development of the GEM- 
3Q* hardware. These algorithms include EMI models of ground response including the effects of a 
non-flat surface/air interface, methods to calibrate EMI instruments based on the magnitude of the sec- 
ondary field from canonical objects, models of the broadband EMI response of ellipsoids, extensions of 
the Standardized Excitations Approach (SEA), and improvements to the Normalized Surface Magnetic 
Source (NSMS) model. Details of these models and algorithms can be found in Sec. 4. 

The GEM-3%" is now a fully functioning EMI instrument capable of vector sensing of magnetic 
anomalies while being well located within a limited range. For precise, queued interrogation of anoma- 
lies, the GEM-3%* provides diverse, accurate, frequency domain data of the secondary EMI field suit- 
able for inversion and discrimination with high fidelity, rigorous models. 
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1 Objectives 


e Develop innovative vector (multi-axis) handheld UWB (ultra-wideband) electromagnetic in- 
duction sensor in the frequency domain, with precise 3-D positioning, for close interrogation 
of anomalies, i.e. for a new higher level of UXO discrimination in the vicinity of a note- 
worthy magnetic responses identified either by the new instrument or most likely by other 
surveying. 


e Develop clutter-tolerant signal processing for UXO discrimination using the data provided 
by the new sensor, based on new, high-fidelity, physically complete forward modeling (the 
Standardized Excitation Approach), rigorous instrument characterization, and on new pro- 
cessing techniques (magnetic field pattern matching, general inversion for equivalent mag- 
netic charge capacity, optimization by Differential Evolution, and classification with Support 
Vector Machines). 


This project builds on SERDP SEED work (project #1353) in response to the needs expressed 
in UXSON-05-03, by developing innovative sensors and signal processing. The aim was and is 
to achieve a new level of UXO survey target characterization and false alarm rejection. The work 
will allow the Army and other DoD components to carry out mandated UXO cleanup more safely, 
effectively, and efficiently. The new fast processing and the portable nature of the products might 
also allow application in tactical situations and for active range maintenance, as well as for cleanup 
away from areas of conflict. 

The positioning and measurement systems here are intended to link to dual mode navigation as 
may be developed under other SONs (Statement of Need). That is, these devices would provide 
the high density, fine mode, while positioning reference points are linked to locations identified in 
the coarse mode. More generally, the survey and processing systems contemplated are designed 
to be most effective under some degree of cuing. The default assumption 15 a high level of cuing, 
in which prior information indicates the approximate locale around an individual anomaly, on the 
order of meters. In both the small and larger scale areal coverage, surveying with a low level of 
cuing is also possible. That is, the cuing would simply consist of identification of some specific area 
some 10’s to +100 m across, within which the new system would accomplish both coarse detection 
and fine interrogation. Reception of the horizontal magnetic field components would help indicate 
the direction of offset from a target location; and the user would then perform arbitrarily detailed 
sweeps around the indicated anomaly. Whether the system to be developed is employed in both 
coarse and fine mode of surveying or only in the latter, the primary focus of the work is on UXO 
discrimination during close interrogation of anomalies. 

This work is designed in part to address realistic conditions: terrain with lumps, bumps, moguls, 
trenches, embankments... and possibly sparse vegetation. High density UWB EMI data will be 
obtained from a new man-portable, vector, FD (frequency domain) instrument with precise posi- 
tioning (sub-cm local precision). The new system is designated below as the GEM-3Dplus (GEM- 
3Y*), meaning a fully integrated GEM-3D plus positioning. The instrument will be “handheld” 
in the sense that a reasonably hardy surveyor will be able to move the sensor head about freely 
in sweeps, tilts, and elevations, for perhaps a couple of hours, without relief. Ergonomic consid- 
erations are not given top priority here, as the result will be primarily a research instrument. The 
new data features provided by the instrument will be exploited by a new generation of original 
processing approaches based on the SEA and NSMS, as well as on other approaches described 
below. The potential benefits envisioned in some of the proposed processing approaches depend 
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on the new level of data content. Also, conversely, the precision and diversity of data contemplated 
only make sense if a new generation of processing can take advantage of it. Thus we will develope 
the instrumentation, survey techniques, and processing together. 

Relative to the system produced during the antecedent SEED work, the new GEM-39%* will 
incorporate a number of desirable developments as listed below: 


1. New small (~ 5 x 8 cm), stackable DSP’s (Digital Signal Processor), each 1 cm in thickness, 
so total stacked height for 3 DSP’s © approx 3 cm. The DPS’s synchronization, power, 
and communication will all be on a bus line (no wiring) with only one communication line, 
and a single RS485 connection at 1MBit/sec to an RS485/USB2 converter as the link to a 
computer; 


2. Each DSP will provide an additional RS232 port for GPS (Global Positioning System) or any 
other external device such as the laser positioning system; 


3. The frequency range will be increased from 48KHz to 96kHz, with low frequency limit of 
30 Hz; 


4. Larger diameter sensor head (54 cm) for more uniform primary field and greater depth of 
view, without significant loss of portability; 


5. The 3 DSP’s will be configured to increase the system performance (convolution is done in 
parallel for all three channels); 


6. The transmitter will be a separate module connected to the DSP via RS422 digital audio 
communication; 


7. New higher power transmitter with four times the voltage of the older system. 


8. Standardized interface and software for the positioning system connection, allowing EM 
(ElectroMagnetic) data synchronization with timing pulses from GPS, inertial, or laser posi- 
tioning; 


9. Integration of positioning system of choice with the EM system so that EM and position data 
files are merged; 


10. The device will be easily reproducible. 


The minimum CPU required for system control will be on board, while wireless links to other 
computers will serve for data storage and processing. This will allow quasi-real time application 
of the processing, including possibly parallelized processing, which can then guide the testing in 
a flexible way. In the least, it will streamline design by separating the control and the storage/pro- 
cessing functions. 


2 Background 


Detection and especially discrimination of buried unexploded ordnance (UXO) from omni-present 
metallic clutter is a persistent, expensive, and pressing problem. While it is not yet clear whether or 
which other technologies might be best used in conjunction with electromagnetic induction (EMI) 
for UXO discrimination, it is clear that EMI is currently a front runner in the development of 
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new, more effective approaches. Virtually all EMI sensors transmit a ’primary” field and receive 
signals ("secondary”’ field) using wire coil antennas. Frequency domain (FD) EMI systems use 
such coils to transmit (Tx) a “primary” magnetic field with selection of frequencies over a chosen 
band - possibly only a single frequency - simultaneously receiving (Rx) and recording responses 
to those frequencies (the ’secondary” field). Thus a central problem in frequency domain sensors 
is the isolation of the receiving coils from the primary field, the latter being very large relative 
to the secondary field, particularly when the receivers are near the transmitters, as contemplated 
here. Physical separation of the Tx and Rx coils is a possible approach. However this typically 
degrades spatial resolution and handy manipulation of the instrument. Therefore, remaining with 
co-located Tx and Rx coils, we implement a “bucking” coil, which transmits a field that opposes 
the primary field at the location of the receiving coil. By canceling out the primary field at the 
receiver as much as possible, one is left (ideally) with only the secondary field from the target. 
Other approaches are possible for finding the secondary field amidst the much stronger primary 
field, e.g. using differencing schemes between stacked or symmetrically distributed coils, but the 
fundamental problem is the same. 

The “bucking” discussion above is particularly pertinent when one attempts 3-D vector defi- 
nition of the secondary field. The geometrical arrangement of the main and the bucking coils is 
crucial to successful suppression of the primary field at the receiver. However, to our knowledge 
no FD EMI instruments have been produced to date that are capable of suppressing the primary 
field in one receiver orientation while also suppressing it in other directions as well. Most imag- 
inable arrangements of two bucking coils that work well for one receiver orientation would in fact 
exacerbate the problem for other receiver orientations. This observation would appear to favor 
the development of time domain (TD) systems for obtaining vector EMI response. Most TD sys- 
tems operate by transmitting a steady signal, saturating a metal object of interest, then shutting the 
transmitter off. A receiver then records the secondary field from the object as that field decays in 
response to the sudden shut off. The receiver only operates while the transmitter is off, thereby 
dodging the primary field. While this bypasses the bucking problem, a number of things motivate 
us to proceed in the FD: FD systems claim the advantage in having superior control of selection 
and power in the frequency content of data produced, and thus in the equivalent time domain sig- 
nal, which can readily be obtained from the FD data. In practice, FD systems are less band limited 
than TD systems and can therefore offer a greater equivalent time range of response than the actual 
TD systems. Consider Fig. 2.0.1 below. 

While there has been progress in ’very early time” EMI devices, in general TD devices are 
unable to shut off the transmitter and damp the effects of the transition quickly enough to get very 
early time data. Information in such data is the equivalent of high frequency FD data. The latter 
is of strong interest because it provides information on asymptotic limits of scattering behavior: 
The induced currents penetrate the target negligibly and the secondary fields they produce do 
not depend on the type of metal encountered. This may offer hope of avoiding very substantial 
complications in UXO discrimination, given that many if not most UXO are composites of different 
metals. At the other end of the spectrum, TD signals naturally fade as the target’s responses wind 
down. This means that very late time data is unlikely to be detectable above the noise. Thus the 
equivalent of low frequency asymptotic limits will be poorly defined or absent. While FD systems 
must also struggle with some challenges in the equivalent very low frequency range, on the whole 
it is much more feasible to get data from them in that part of the EMI band. Many UXO produce 
significant response patterns in the 10’s of Hz. Low frequency asymptotic limits and fundamental 
signal features containing basic target shape information appear in many cases only well below 30 


Benjamin Barrowes, Kevin O’ Neill -3- - ENGINEERING 
Ἐς: B44) AT DARTMOUTH 


MM-1537, GEM-3 2): -- Final Report 2 BACKGROUND 
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Figure 2.0.1: General representation of received time domain signal (R) vs log of time. 


Hz [SERDP Project 1122, Final Report]. To get at the TD equivalent of these frequencies down to 
1 Hz using TD instruments would require recording out to the order of tenths of a second. Highly 
unlikely. At least at the research level, FD EMI measurements have been made successfully over 
a continuous UWB ranging from a few Hz up to about 300 kHz, roughly the TD equivalent of 
spanning times from a micro-seconds to tenths of a second. Recent work by Geophex on TD 
systems in a GEM configuration may be an exception to the limitations cited above, and we look 
forward to encountering that new device. 

Much of the motivation for this project derives from a desire to develop further a new generation 
of EMI forward modeling; to apply it in innovative processing; and to develop instrumentation that 
will provide the data most likely to support this new processing to best advantage. The new forward 
modeling will generally be based on the Standardized Excitations Approach (SEA), though all 
processing will not be restricted to SEA implementations. The theory and detailed development 
of SEA have been presented in various fora [1-12]. Here we only present a general sketch of the 
underlying ideas so that the direction taken in some of the proposed work may be understood. 
The essence of the method lies in its ability to provide physically complete forward modeling - 
including all target heterogeneity, near and far field, and internal interaction effects - while being 
fast enough for inversion and classification computations. 

During antecedent SERDP SEED work, the configuration of the existing Geophex GEM-3 [4] 
was used as a starting point for the new vector device. The essential difference was that the new in- 
strument would receive and record separately three orthogonal magnetic field components instead 
of one. This was accomplished by adding two receiver coils that are the same size as the original 
single coil but which are perpendicular to it and to each other. These additional coils escape be- 
ing swamped by the primary field in part because they are near the null region engineered for the 
original in-plane receiver loop. For the most part, however, they are protected by their orientation. 
The plane of each of the new receiver coils is completely parallel to the transmitted field lines from 
the perpendicular transmitter coils (Fig. 2.0.2). Because the receiver coils respond only to rate of 
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change of magnetic flux that passes through them, in principle the additional coils do not register 
the primary field. The bench level instrument that was constructed and tested has an expanded 
bandwidth relative to established versions of the GEM-3, ranging from the same lower limit of 30 
Hz up to at least 98 kHz, with a maximum of about 15 frequencies recordable in each channel, i.e. 
for each orthogonal signal component (using slightly fewer is recommended). 


receiver 


Transmitter 


Figure 2.0.2: An EMI sensor field impinges on a UXO, which tends to respond most strongly along its axis, 
producing a response that strikes the receiver at an arbitrary angle relative to the vertical. 


The bench version of the instrument that was produced appears to work well by both quan- 
titative and qualitative criteria. Proper (and illuminating) symmetries and asymmetries appear in 
the different received field components as the sensor head is swept over UXO’s and other objects. 
Measured data for machined spheroids of different materials and shapes compare very well to cor- 
responding analytical solutions that account for the details of the GEM’s primary field. The laser 
positioning system that was implemented is said (by the supplier) to be capable of sub-millimeter 
accuracy, and was configured to provide all tilt angles in a version of the instrument dubbed the 
GEM-3DL. In various benign tests the integration of the system performed well. 

Regarding the EMI measurement system, the task that remains is proper design and imple- 
mentation so that it is melded into an efficient and coherent unit, including direct integration with 
the positioning system via a conventionally recognizable timing pulses. Basically, for the sake of 
showing that the basic concept of FD vector EM measurement could work, the SEED work ef- 
fectively patched together in parallel three traditional GEM-3’s, one associated with each receiver 
loop. Given what was learned during the SEED exercise and also the recent advancement in com- 
ponents at affordable prices (e.g. digital signal processors — DSP’s), we are now positioned to 
construct and integrate the instrument as it should be done. It has also become clear that the easiest 
and most reliable way to perform the EM-positioning integration is to key the EM measurement 
timing to a pulse from the positioning system, as opposed to vice versa. 

The GEM-3D produced in the SEED effort used a laser positioning system: the ArcSecond 
laser positioning system. The laser positioning system provided accuracies on the order of a few 
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mm of the instrument head. The positioning data provided by this system allows free movement of 
the sensor head at any elevation and tilt angle thus leading to a rich and diverse data set. Data from 
multiple elevations has helped in distinguishing target signals from widespread clutter responses; 
and in the general inversion/optimization calculations we have done, it significantly reduced or 
eliminated ill-conditioning. Also, horizontally distributed observation points allow “multi-look” 
data, i.e. views that stimulate and also register response from all sides of an object to be studied. 
All this requires very flexible positioning, 1.6. both antenna movement and tracking thereof. In fact, 
without such positioning the new (and older GEM-3) sensor will not be useful for the next gen- 
eration of processing towards which we strive. However significant design and integration issues 
remain. Both inadequate calibration of the laser system and especially timing disparities between 
both the different laser receivers and between the laser and the EM system precluded continuous 
motion of the sensor in recent field tests at the Aberdeen Proving Ground UXO Standardized Test 
Site. This makes for impractically slow surveying, even for testing purposes. 

Initially, we planned to use the ArcSecond system with the new GEM-3Q%* in order to achieve 
the same level of positioning accuracy. However, early on in the project we reconsidered this design 
decision in light of recent developments in the area of “beacon” positioning systems (see Sec. 3.3) 
and because of factors mentioned above regarding the GEM-3Y*. Beacon positioning systems 
rely on the primary field of the instrument itself to provide an accurate position of the sensor head. 
Using external receivers to detect and sample the primary field of the GEM-39*, the location 
of the primary field (assumed dipolar) source can be inferred. The design, implementation, and 
resulting accuracy of this system is detailed below. This positioning system, in theory, is simpler 
to operate and maintain than the ArcSecond system. For example, the calibration of the system is 
far simpler for the beacon system, and the setup time is also greatly reduced. On the other hand, 
the range of this system is less than the laser positioning system. 

In simulations and lab measurements, we have seen distinct benefits from vector data in the 
task of defining unambiguously the response of UXO’s and of other objects to the basic modes of 
excitation. Vector or “3-D” data means that three orthogonal components of the scattered magnetic 
field are measured at the receiver location, thereby completely defining the field there that consti- 
tutes the object’s response. The idealized representation in Fig. 2.0.2 shows how a UXO response 
field will generally strike the receiver at some arbitrary angle relative to the vertical, i.e. relative to 
the direction in which the established sensors record the signal. 

The rest of the report proceeds as follows. Sections 3 and 4 detail the materials and meth- 
ods associated with the design and implementation of the hardware of the GEM-39* as well as 
the software analysis algorithms developed in the course of this project. Section 4 includes dis- 
cussions on some characteristics of the GEM-3Y* (Sec. 4.1), the soil response of the instrument 
(Sec. 4.2), absolute calibration of the instrument (Sec. 4.4), the Generalized Standardized Excita- 
tions Approach (GSEA, Sec. 4.5), the NSMS and HAP methods (Sec. 4.6), and the response of 
the GEM-39* to ellipsoidal targets (Sec. 4.7). These sections are followed by blind test results 
both from lab data and from test plot data at CRREL (Sec. 5). The Appendices contain supporting 
material such as operating procedures and a publication list follow the concluding section (Sec. 6). 
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3 Materials and Methods: Hardware Development of GEM-39* Instru- 
ment 

3.1 GEM-39* Design 

3.1. GEM-39~ Specifications 


The primary purpose of this project is to advance the existing GEM-3 EMI sensor to incorporate 
multiple axis viewing (’vector”’) capability that is expected to enhance the UXO detection and 
discrimination. Relative to the system produced during previous CRREL project, the new GEM- 
3Z* will incorporate a number of advancements as listed below: 


e Larger diameter sensor head (54 cm) for more uniform primary field and greater depth of 
view, without significant loss of portability 


e The frequency range will be increased from 48 kHz to 96 kHz, with low frequency limit of 
30 Hz 


e Use of new small stackable DSPs, each about 1 cm in thickness. Synchronization, power, 
and communication will all be on a bus line. An RS485 connection at 1Mb/sec to an 
RS485/USB2 converter will serve as the link to a computer 


e Each DSP with an additional RS232 port for GPS or any other external device such as the 
laser positioning system 


e DSPs configured to increase the system performance (convolution is done in parallel for all 
four channels) 


e A separate transmitter module connected to the DSP via RS422 digital audio communication 
e New higher power transmitter with increased drive the voltage 


e Standardized interface and software for the positioning system connection, allowing EM data 
synchronization with timing pulses from GPS, inertial, or laser positioning 


e Integration of positioning system of choice with the EM system so that EM and position data 
files are merged 


e Easily reproducible. 


3.1.b Early Progress 


This section describes our early design efforts for the GEM-3Y* sensing head, transmitter and 
receiver coils, and structural assemblage of the coils and mechanical support structures. The GEM- 
3Y* has two additional orthogonal vertical coils in addition to the existing GEM-3 geometry. 
Figure 3.1.5 shows the overall appearance of the GEM-39Y* sensing head. The design provides 
the following features: 


e Sensing head: 1/2” thick double-skin fiberglass foam core synthetic board with an overall 
diameter: 55 cm; 


e Outer transmitter: 54 cm in diameter; 
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e Receiver coils: one circular horizontal coil in the magnetic cavity plus two rectangular ver- 
tical coils. The vertically offset rectangular shapes reduce the mandatory sensor ground 
clearance 


e Maximum visibility of the target area by embedding many see-thru holes in the disk 
e Minimum sensor head weight; 


e The electronic console and battery to be mounted at the end of the handle bar, partly to act 
as a counter weight. 


All three receiver coils, one circular and two rectangular, are designed to have an identical 
area-turn product. The two rectangular coils, as shown in Fig. 3.1.5, are symmetrically centered 
across the sensing disk. The main advantage for centering is that all three receiver coils are mathe- 
matically co-located (this was later changed, see below), which simplifies locating and interpreting 
targets. Otherwise, these coils may move arbitrarily upward. 

Figure 3.1.1 shows the design sheet for this GEM-39*. Basic design parameters include the 
diameters, turns, and wire sizes for Tx1 (outer transmitter), Tx2 (inner bucking transmitter), and Rx 
(the horizontal receiver coil). The two rectangular vertical coils (RxV1 & RxV2 -not shown in this 
design sheet) are made such that their area-turn product is identical to that of the horizontal coil. 
This sheet is based on the 12VDC supply, which is a mere scale factor for the transmitter moment. 
Higher supply voltages are used in the final design. The design shows a healthy transmitter moment 
of about 35 A-m2 at 750 Hz, which compares favorably with the previous 40-cm diameter version 
having about 10 A-m2 also at 750 Hz. The dipole loss due to Tx2 is only 5.37%. The sheet also 
shows the I- and Q-responses of the Rx in ppm to a known Q-coil between 30 Hz—96 kHz. The 
specs and location of the Q-coil are also indicated in the sheet. 

Figures 3.1.2—3.1.6 show some preliminary drawings of the GEM-39*. Later designs posi- 
tioned the transverse receiver coils to be more flush with the bottom of the transmitter coil (see 
Fig. 3.2.6). 


3.2 GEM-3%* Continuing Development 
3.2.a Cord Sensitivity Issue 


During the initial stages of testing the GEM-3Y* at CRREL, we noticed an extreme sensitivity in 
the real part of the transverse components associated with any movement of the instrument. We 
eventually traced this sensitivity problem more specifically to issues of cord movement. Whenever 
the connecting cord from the sensor head to the electronics box was moved in the slightest amount, 
the transverse components of the received magnetic field were corrupted by a tens if not hundreds 
of ppm. Together with the Geophex engineers, the problem was traced back and solved so that this 
sensitivity issue was resolved. 


3.2.b Head Deformation Issue 


Prototype grade instruments often have several unforeseen issues arise during their development. 
After the sensitivity associated with cord movement was resolved, we noticed another type of sen- 
sitivity issue whenever we physically touched or moved the instrument. As it turned out, there was 
a further sensitivity associated with deformation of the relative position of the transverse receiver 
coils to each other. We decided to investigate whether it was the deformation of the entire head or 
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Handle Bar 


Figure 3.1.1: GEM-3Y* basic geometrical design. Table showing coil parameters and frequency character- 
istics omitted: proprietary. 
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Figure 3.1.2: ΘΕΜ-3 929 preliminary drawing 1. 
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Figure 3.1.3: ΘΕΜ-3 929) preliminary drawing 2. 


Benjamin Barrowes, Kevin O’ Neill -l1- 


THAYER SCHOOL OF 


ENGINEERING 
AT DARTMOUTH 


MM-1537, GEM-3 9*-— Final Report 3. MATERIALS AND METHODS: GEM-39* HARDWARE 
3.2 GEM-3%* Continuing Development 


GEM-3D+ RX coil form 
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Figure 3.1.4: GEM-3Y* preliminary drawing 3. 
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Figure 3.1.5: GEM-39* preliminary drawing 4. Final design for receiver coils is shown in Fig. 3.2.6. 
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Figure 3.1.6: GEM-3Y* preliminary drawing 5. 
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Figure 3.2.1: Pictures of the GEM-3Y* in development at an intermediate stage before finalization of the 
sensor head. 
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just of the receiver coils that was causing this broadband shift. So, we squeezed the two transverse 
receiver coils toward each other, without touching the main head, and sure enough we got only in 
phase broadband shifts up or down depending on which direction we squeezed in. 

Any slight misalignments from perfectly perpendicular to the primary field should be consid- 
ered background and subtracted if it was constant. The problem is that these transverse receiver 
coils are not as rigidly fixed in their geometry with respect to the primary field as the Z receiver coil 
is. The transverse receiver coils are only embedded in a 1 cm square cross sectional scaffolding 
that is much less robust in terms of torsional and linear deformation. Similarly when I squeeze the 
receiver coils toward each other, even if I deflected them only by a few microns, it has this large 
effect of up to 100 or 200 ppm in the data. The fact that these deformations have a broadband, con- 
stant (with respect to frequency) inphase shift leads me to suspect it is the primary field response I 
am seeing. 

To lessen this effect in the short-term we simply created some handholds on the supporting 
plexiglass structure we now have on the GEM-3%* so that we can pick it up from the center bot- 
tom of the instrument without putting any torque on the main head (see Fig. 3.2.2). We ended up 
stiffening the head with plexiglass as shown in Fig. 3.2.2. This helped the problem considerably, 
but it neither looked good nor was it very portable due to the increased weight. The engineers at 


ae 


Figure 3.2.2: Plexiglass gussets glued to GEM-3QY* sensor head. 


Geophex redesigned the GEM-39* sensor head using lightweight foam core epoxy coated stiff- 
ening members. A current picture of the GEM-39~ including these and members is in Fig. 3.2.3. 


3.2.b.(1) Final Head Design Figure 3.2.6 shows the final design of the GEM-39" sensing head. 
The only difference from the design above is that the two vertical coils (Hx and Hy) are raised 
by 3.5 cm, mainly to reduce the ground clearance distance. This was done at the request of the 
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Figure 3.2.3: Foam-core gussets epoxied to GEM-39* sensor head. 


CRREL project manager, Dr. Ben Barrowes. Figures 3.2.4 and 3.2.5 show the sensor head during 
the construction and initial wiring stages. 


3.2.b.(2) Sensor Electronics Figure 3.1.1 is the final design sheet for this GEM-39*. The main 
change from the previous report includes the transmitter wire gauge that changed from 14 strands 
of AWG#30 to a Litz wire of 110 strands of AWG#30, substantially reducing its resistance (from 
0.251 to 0.032) but adding about 0.45 kg (1 pound) of copper weight. The combined transmitter 
(TX1 + TX2) shows a measured inductance of 53H (close to the predicted value) and a self- 
resonance frequency at 1.44 MHz (highlighted yellow in the upper right corner (of proprietary 
table)), well above the intended operating bandwidth. 

The wire gauge of the three receiver coils has also been changed from AWG#30 to AWG#26 
to reduce the resistance. The three coils, designed to have similar area-turn products, have the fol- 
lowing measured electrical parameters (Table 1). Small variations are within acceptable tolerance 
and will be equalized through calibration procedures, once the sensor has been completed. 

We assembled the sensor electronics beyond the sensor head in year 2 of the project. Basically, 
there is a single electronic console mounted at the end of the handlebar, partly acting as a coun- 
terweight for handheld operation. The console, including a rechargeable battery, communicates 
wirelessly thru Bluetooth with either a PDA or a laptop that will display data and graphics during 
a survey as well as archive the raw data. The Bluetooth also communicates with the GEM Beacon 
Navigator (fabricated in Year 2) to receive wirelessly the sensor location data. 

Figure 3.2.7 shows the overall system functional block diagram for the GEM-39*. The re- 
ceiver diagram, Fig. 3.2.8, utilizes a custom-designed DSP board that, among others, performs 
discrete sine- and cosine-convolution operations at each transmitted frequency to obtain the ppm 
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Figure 3.2.4: GEM-39* sensing head and coil frames 


Figure 3.2.5: Beginning stage of wiring. All three receiver coils are now potted, electrically shielded (silver 
coating), and epoxy-painted for protection. Notice a small preamp printed circuit board (PCB) for each 
receiver coil taped to the round z-axis receiver coil. 


Benjamin Barrowes, Kevin O’ Neill -18- - ENGINEERING 
Bs 6/4) AT DARTMOUTH 


MM-1537, GEM-3Y*- Final Report 3. MATERIALS AND METHODS: GEM-39*+ HARDWARE 
3.2 GEM-39%* Continuing Development 


Handle Bar 


Figure 3.2.6: Final drawing of the CRREL GEM-3%* sensing head. 


mean coil offset from #of measured measured _ self-resonance 
coil dimensions TX plane turns resistance inductance frequency 
Hz 7.25cm dia circle 0 100 6.20 2.78mH 254kHz 
Hx | 10.3cm x 14.6cm rect +3.5cm 110 7.22 3.07mH 242kHz 
Hy | 10.4cmx14.7emrect +3.5cm 106 7.0Q 2.92mH 242kHz 


Table 1: Receiver coils characteristics 
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responses of the inphase and quadrature components. Figure 3.2.9 shows the transmitter functional 
diagram. 


Bit Stream 
and Clock 
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Computer a Probability 
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Figure 3.2.8: Custom-designed DSP module for the GEM-39*. 


3.3 GEM-32%" as Beacon Positioning System 


The GEM-3%* acquires high quality data from its three receiver coils. However, for a hand- 
held instrument such as this, accurate positioning of the sensor head including all three Euler 
angles is very important. One reason why position is so important is that the depth to which the 
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Power Supply 


Digital 
Audio 
RX 


Bridge Drivers 


Figure 3.2.9: Transmitter block diagram with red lines showing power flow and orange lines data flow. 


GEM-3%* can see” is limited to approximately 70 cm for realistic size targets. At the shallow 
depths, a positioning error on the order of 2 cm can corrupt data and subsequent inversion schemes. 
Accuracy of all three orientation angles is important for translating local vector components to an 
invariant global system during rotations. As a consequence, we expended significant efforts at 
developing and integrating the beacon positioning system into the GEM-39°. 

Most currently used positioning systems rely on locating something attached to the sensor, 
for example the robotic total station, GPS receivers, and laser positioning systems. The beacon 
positioning system we have integrated into the GEM-39* works on a different principle than 
these types of positioning systems. It makes the assumption that the primary field transmitted 
from the sensor head is dipolar in shape. Under this assumption, two triaxial receiver coils in 
a fixed configuration should theoretically allow the determination of the location and two of the 
three Euler angles associated with the sensor head. Because the dipole field is axially symmetric, 
the third Euler angle which describes the rotation of the sensor head around its axis of symmetry 
cannot be deduced simply from the dipolar primary field. In fact, for the usual GEM 3 sensor head, 
this last Euler angle is unnecessary. However, because the GEM-3Y* has the additional transverse 
receiver coils whose responses depend on this rotation about the axis of symmetry, some method 
to find this last Euler angle was needed. We decided to install an electronic compass onto the main 
sensor head. This compass when combined with the data from the beacon positioning system, will 
provide the full set of Euler angles necessary to uniquely determine the location and orientation of 
the sensor head at all times. 

The beacon positioning system that Geophex demonstrated to us in 2007 was for a different 
Geophex instrument. This other beacon positioning system had a constant electronic gain built 
into the system itself. This gain is an important factor in determining the sensitivity and accuracy 
of the final positioning system. Too small of a gain leads to a short range; too large of a gain will 
lead to saturation when the sensor is near the beacon’s triaxial receivers. For the GEM-39", the 
decision was made to have a variable gain in order to have both a large range from the receiver 
coils and accuracy when the sensor head is in close proximity to the beacon receivers. This added 
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Figure 3.2.10: Detail of the transverse receiver coils on the ΘΕΜ-3 2).. 
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Figure 3.3.1: Concept of GEM-39* positioning system. The (assumed) dipolar shape of the GEM-39* 
primary field is recorded by two triaxial receivers held in a fixed relationship. 


complexity to the positioning system delayed its final integration with the GEM-39~. 

Conceptually, the GEM-3Y* beacon positioning system operates as depicted in Figs. 3.3.1 
and 3.3.2. The data from the beacon system is fully incorporated into the GEM-39Y* data stream. 
One important note about this positioning system is that the EMI primary field from the GEM-39* 
penetrates trees and mounds as easily as it does the soil, so this system can operate in challenging 
environments with non-metallic obstacles. 

A schematic of the coordinate systems involved in the GEM-39* beacon system is shown in 
Fig. 3.3.3. 


3.3.a Straight Line Tests 


As an initial test of the beacon positioning system, we set up a 2x10 on the ground at the Geophex 
facility and moved the GEM-3 along straight lines next to the boards at different distances from 
the beam. Figure 3.3.4 shows the setup used in this test, while Figs. 3.3.5—3.3.8 shows the relative 
accuracy of the system at near, intermediate and far distances from the beam. The processing 
system infers the x, y, and z position of the sensor head as the instrument is moved around the 
plank. The merged data, EMI and positional, are then recorded simultaneously. If desired, the 
position of the sensor head in real time can be displayed (see Fig. 3.3.9). This display could be 
useful for real time detection of anomalies as both the EMI data and the sensor head position are 
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Figure 3.3.2: Beam holding the two triaxial receivers of the GEM-3Y* positioning system. 


Beacon positioning system 


Figure 3.3.3: Coordinate systems for the target, beam, and GEM-397. 
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Figure 3.3.4: Beacon positioning system and plank setup for initial demonstration at Geophex facility. 
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Figure 3.3.5: Beacon positioning system data along a board for the entire length. 
overlaid. 
The details of the beacon positioning system involve the nonlinear inversion of a dipolar source 


based on 6 measurement points. The instrument and the two triaxial receivers are shown concep- 
tually in Fig. 3.3.10. An oscillating dipole magnetic field is generated by the EMI transmitter. The 
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Figure 3.3.6: Beacon positioning system data along a board. Near range. 
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Figure 3.3.7: Beacon positioning system data along a board. Intermediate range. 
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Figure 3.3.8: Beacon positioning system data along a board. Large range. 


Figure 3.3.9: Beacon positioning system real time feedback of EMI data (not shown) and position. 
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Figure 3.3.10: GEM-39* beacon positioning system principles. 
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where R = \/x*+y?+2? dipole-field point distance, Uo is the permeability of free space, the 
subscripts refer to field component/source dipole component, and the measurement (field) point is 
at the origin. For the locator system comprising two sets of coils at two field positions (xa,0,0) 
and (xb,0,0), two sets of equations as above are formed by replacing x with (xa,xb). We re- 
parametrize these equations in terms of a total dipole strength and the sine of the tilt angels in 
place of m,.y,, A non-linear iterative inversion algorithm is performed by a computer to solve For 
the x,y,z and tilts; at initialization we also solve for the dipole strength, then freeze it for stability. 

The accuracy of the beacon positioning system was determined in tests at the Geophex facility 
in 2007. The complete setup can be seen in Fig. 3.3.11. The sensor head was then moved on a 
precise grid while the beacon positioning system was recording data. The calculated position was 
then compared to the actual position on the grid with results shown in Fig. 3.3.12. The horizontal 
projection of locator measured relative position (blue diamonds) vs. template grid (red crosses). 
There are 50 samples at each position. The 3-dimensional RMS error is 0.77 cm. Note, one grid 
point is missing due to data acquisition dropout. Figure 3.3.13 repeats the test except the nearest 
edge of the gridboard was 1.5m away from the sensor head. Relative measured position (blue) 
plotted on the ideal grid (red) for the template at 1.5 m closest edge from the locater beam. There 
are 50 samples superposed at each point. The 3-D rms error is 0.92 cm. 

More information on the processing method used in analyzing dynamic data from the GEM- 
342) } beacon positioning system can be found in Sec. 4.1.c.(3). A more detailed explanation on the 
construction of the beacon hardware can be found in Appendix C.2. 
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Figure 3.3.11: Complete setup of the beacon positioning system tests at Geophex. 


3.4 GEM-3%* Complete 


Figure 3.4.1 show the completed GEM-39* before the stiffening was added (see below). Also see 
Fig. 4.1.1 for picture of the completed GEM-3Y* instrument head. 

The GEM-39%* instrument is now completed, including the receiver beam for the beacon po- 
sitioning system. Software which runs both the beacon software and the GEM-39" itself has 
been provided to CRREL by Geophex. Setup and operation is straightforward and is detailed in 
Appendix C.1. We have acquired both lab and test plot data and run blind tests of these data (see 
Sec. 5). As with most prototype instruments, the GEM-39* is not a production ready instru- 
ment. There remain some issues with software stability and hardware communication (between 
the GEM-39* and the beam). These issues cause software crashes and resets that restrict larger 
deployment. The head deformation issue described above is still relevant when the GEM-39* is 
held in a dynamic data acquisition mode. Essentially, when the GEM-39* is tilted beyond a nomi- 
nal degree (beyond ~ 10°), the primary field can leak into the transverse receivers and contaminate 
the data. 

Frequency domain instruments such as the GEM-39%* are desirable for several reasons cited 
above (Sec. 2). Lessons learned from the GEM-39%* can assist future researchers in developing 
more sensitive and robust FD instruments and beacon positioning systems that can potentially 
resolve these issues. 
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Figure 3.3.12: Results for beacon positioning system test #1. The nearest edge of the gridboard is 1m from 
the beam. The horizontal projection of locator measured relative position (blue diamonds) vs. template grid 
(red crosses). There are 50 samples at each position. The 3-dimensional rms error is 0.77 cm. Note, one 
grid point is missing due to data acquisition dropout. 
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Figure 3.3.13: Results for beacon positioning system test #2. The nearest edge of the gridboard is 1.5m 
from the beam. Relative measured position (blue) plotted on the ideal grid (red) for the template at 1.5 m 
closest edge from the locater beam. There are 50 samples superposed at each point. The 3-D rms error is 
0.92 cm. 
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Figure 3.4.1: Final GEM-3Y* instrument. 


4 Materials and Methods: Software Modeling and Analysis of GEM-39~* 
Data 


During the course of MM1537, we developed models and algorithms which aid in helping to 
analyze not only the data from the GEM-39", but potentially from any other EMI instrument. This 
section will proceed as follows. Some general characteristics of the GEM-3Y~* and the beacon 
positioning system are presented as Sec. 4.1. Also in that section is a preliminary description of the 
GSEA (Generalized Standardized Excitations Approach), while a more detailed presentation is in 
Sec. 4.5. In Sec. 4.2, we first discuss the response of two different soils to the GEM-3Y* acquired 
at a local gravel pit. We also discuss the form of the permeability of most soils at EMI frequencies 
and discuss implications for both FD and TD instruments. We also model soil properties such 
as surface roughness and discuss the estimation of soil susceptibility in Sec. 4.3. Section 4.4 
documents our efforts to find an absolute scale factor of GEM ppm converting it to magnetic 
field H values. This can be useful for local soil characterization and background subtraction, in 
particular. 

We have adapted our normalized surface magnetic source (NSMS) and GSEA models to the 
GEM-39* geometry (see Sec. 4.6 and Sec. 4.5). These ultra-high fidelity, physically complete 
models facilitate discrimination processing by a) allowing detailed inspection of expected signals 
under any contemplated circumstances; b) enabling signal pattern matching to establish the pres- 
ence of known or likely target types; and c) providing fully realistic synthetic input data on a large 
scale for training statistical learning machines. In addition, we also developed a model for EMI 
responses by elliptical metallic bodies, applicable at higher frequencies by design, but also match- 
ing fairly well with data at low frequencies if they are highly permeable (see Sec. 4.7). This is the 
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Figure 4.1.1: GEM-39* sensor and measurement grid 


only non-numerical solution available for objects with three dimensional geometry. We have also 
machined 6 ellipsoids of varying aspect ratios to compared our models with measured data [5]. 
These models will allow us to accommodate different size and shape targets whether they are in 
our library for not. 


4.1 A vector handheld frequency-domain sensor for UXO identification 
4.1.a Executive Summary 


This section presents some methods and results related to UXO identification using the GEM-39". 
Our analyses exploit data provided by the sensor in both grid-based and dynamic measurements to 
characterize different objects, including metal spheres and actual UXO. For the data analysis we 
alternate between the dipole model and the more rigorous standardized excitation approach. We 
review some ill-conditioning issues encountered with the latter model and the different approaches 
that we use to overcome them. In applications, the availability of horizontal field components in 
the data allow us to identify UXO vs. non-UXO items while minimizing the nonnegligible effects 
of ground response. 


4.1.b Introduction 


The improvement in sensing methods represented by the GEM-3Y~* must be accompanied by fast, 
physically complete, and accurate forward models that can both reasonably pinpoint the location 
and orientation of an object and provide sufficient information on its electromagnetic properties so 
as to characterize it unambiguously. Here we adapt to the GEM-39* some forward models that 
have shown success with other sensors. Perhaps the most widely used model treats each anomaly 
as a point dipole [6—10]; this procedure, while limited, is fast and provides good starting guesses 
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for more elaborate methods. An example of the latter is the data-derived Standardized Excitation 
Approach (SEA), which incorporates the finite size of the scatterer and the nonuniformity of the 
primary field in the relevant length scale. The approach decomposes the primary field striking the 
target into a set of standardized modes, whose relative importance in the expansion is given by 
the particulars of position and orientation. The object is replaced by rings of magnetic sources 
placed on a surrounding spheroid. The response of the rings to each mode of excitation can be 
described by a Green function, and the full response of the target is a superposition of these with 
the corresponding geometric coefficients as weights. Once the object’s response to each mode is 
established, the response can be calculated for any variation of excitation and receiver location by 
just using an appropriately altered superposition of the modal responses. The method was previ- 
ously adapted to the GEM-3 [11] and to various time-domain instruments [12, 13] using charges as 
the responding sources; a generalization replaces the charges by dipoles and uses Stokes’s theorem 
to streamline the necessary integration over the receiver areas [14, 15]. 

The overall section is organized as follows. After describing the instrument in Section 4.1.c we 
present some studies related to UXO identification that use data produced by the GEM-3* in both 
grid-based and dynamic measurements. We employ both the dipole model and the SEA, which we 
introduce in sections 4.1.d.(1) and 4.1.d.(2) respectively, to characterize different objects, including 
metal spheres and actual UXO. We look at the possibility of using the horizontal components of the 
secondary field to identify UXO vs. non-UXO items while minimizing the nonnegligible effects 
due to the ground. We conclude in Section 4. 1-f. 


4.1.5 The GEM-32~* Sensor 


The GEM-39%* sensor, designed and manufactured by Geophex, Ltd., is a wideband frequency- 
domain instrument that continuously transmits chosen frequencies and outputs FD spectra of re- 
sponse. The user can select up to 15 frequencies of operation for the sensor, ranging between 
30 and 90030 Hz. The measured secondary field at each frequency is given in dimensionless units 
of parts per million (PPM), defined as [16, 17] 


PPM — Σουρ ἐΠΆΕΌΤΙΟ flux through the receiver 66: 10%, (4.1.1) 
Primary magnetic flux through a reference coil 


which serve to divide out the frequency dependence and the value of the transmitted current and 
account for the fact that the actual measured quantities are voltages proportional but not identical 
to magnetic fluxes. 

The sensor head contains two circular coplanar concentric transmitting coils connected in se- 
ries. These transmitters have radii 27 cm and 10.84 cm and are wound 6 and 2 turns respectively, 
making the effective current on the outer coil three times larger [18]. This “transmitter bucked” 
arrangement is designed to produce no primary magnetic flux through a receiver coil, coplanar and 
concentric with the transmitters and of radius 7.25 cm, to which we henceforth refer as the H,-coil 
(these numbers are determined by using an expression equivalent to (4.1.2) below but cast as an 
expansion in Legendre polynomials [16]). The H,-coil can then measure the much smaller flux due 
to the secondary field. 

The other two receiving coils on the GEM-39Y* are 


1. The “Ἢ, coil with dimensions 10.3 cm (in 2) by 14.6 cm (in %), offset 3.5 cm from the plane 
of the transmitter coils; 
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20 


Figure 4.1.2: (Left) The field H?" due to the GEM-39* slightly above the plane of the transmitter coils. The 
central cavity where the primary flux vanishes is clearly visible. (Right) The field H!* due to the GEM-39* 
on the plane of the transmitter coils, this time seen in profile. 


2. The “Hy” coil with dimensions 10.4 cm (in 2) by 14.7 cm (in ὦ), offset 3.5 cm from the plane 
of the transmitter coils. 


The primary flux through these vertical receivers vanishes trivially. 
The primary field due to a circular coil of radius x in cylindrical coordinates is [19, 20] 


Be yee £1 (1-34) x(a) £0), where the modulus ΚΖ = “Px ee 
1 k\ip 2 (~+p)*+z 


(4.1.2) 


where K(k) and E(k) are elliptic functions [21] 

The field H!” slightly above the plane of the transmitter coils is depicted in three dimensions 
and as a function of p in Figure 4.1.2. This figure also shows how the total field arises as the sum 
of the contributions of the two transmitter coils. 


4.1.c.11) Integrating over the areas of the receivers Consider an upright dipole located zg m below 


the origin. With m = m2, r = pf, and r’ = -- χρῶ we can use the general formula 
-R)R/R? - 225 — p? 
ὦ δ eS eee pe. (4.1.3) 
47R3 An (2+ p2)5/2 
Integrating (4.1.3) over the area of the horizontal receiver we find 
a 2. 
m a 
Φ-- 2 Π,ράαρ----------. 4.1.4 
τ zp ρ 2 (25 +.a?)3/? ( ) 


On the other hand, if we substitute the field at the center and multiply by the receiver area, we get 
the “average” flux 


2 3 
ei τ πε (4.1.5) 
4π ζῃ 2 ζῃ 
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We see that the average flux is larger than the integrated counterpart when we are right above the 


target: 
Φ' α12]}2 
== 1 | (4.1.6) 


For fixed a this effect grows as zo decreases; i.e., when the target gets closer to the sensor. 
We can also let vary the horizontal position of the target. We take r’ = χρῶ — zo%, and we obtain 


for the field oe 
2% — p* — 2 
ane i δεῖ. τε, a (4.1.7) 
An (29 + p* +x5 — 2pxocos Φ) 
There is probably a closed expression for the flux, but it is undoubtedly very complicated. An 
intermediate expression is 


a 2 13 2. οὐ 4,2 42 2\ 2,4 
oa m ὧν pdp 5 5 K(k) a ((p *0) 2za(P + %0) 3z) e(), 
8x Jo (p|xol)?/?(p — |xol)? +29) zo ((p — |xo|)* + 56) ae 


where K(k) and E(k) are the complete elliptic integrals of the first and second kinds and the 
modulus 
ἔξ τ -- (4.1.9) 
(ρ + |xol)? + τῷ 
Making x9 = 0 in (4.1.8) we easily recover (4.1.4). If we merely substitute the field at the center 
and multiply by the receiver area, we obtain 


a a i δ πα (4.1.10) 
An (ζῇ +x9)>/2 4 (+26)? 


As part of the Generalized Standardized Excitation Approach to be considered later, it is pos- 
sible to use a shortcut in this integration by employing Stokes’s theorem: 


facd= [(vxa)-rds (4.1.11) 
Cc 5 


Defining a “vector potential” (which in rigor should be proportional to the permeability of the 
medium) by 


mx R 
= —_ 4.1.12 
4π}3 ( ) 
we can easily verify that its curl is given by (4.1.3). For the situation considered in this section we 
have " Ε 
ΠΡῸ in generaland A= ae (4.1.13) 


~ 4n(p2 +23? 4m (a? + 22)3/2 


along the circumference of the receiver, which trivially integrates to yield (4.1.4). 


4.1.c.(2) The other receivers To verify that our integration routines are correct we compare the 
integrated flux to the average one for an upright dipole. In this case the depth is z9 = —35 cm, and 
we have xp = 4.cm and yo = 12 cm. Interesting results: The ΗΠ plots are very similar; the Ης plots 
show the same tendency we saw in the last section—the averaged flux tends to be larger than the 
integrated one; the Hy plots show the opposite tendency. See Figure 4.1.3. The dipole can also be 
tilted, as Figure 4.1.4 shows. 
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Figure 4.1.3: Integrated (blue) and averaged (green) flux profiles on a grid above an upright dipole; from 
left to right we have the z-, x-, and y-components. The red is the difference between the two. 


4.1.c.(3) Beacon positioning One of the features that make the GEM-39* potentially useful for 
dependable UXO identification is its positioning system, which provides sub-centimeter accuracy 
and can be used in treed or rugged terrain. The beacon is also lightweight and unobtrusive, requir- 
ing almost no hardware to be attached to the sensor itself. 

For the data analysis it is necessary to have the tilt information referred to the same global 
coordinate system in which the sensor location is known. This cannot be inferred directly, however, 
since the field of a dipole is azimuthally symmetric: the GEM-39* can be rotated about its z-axis 
and appear identical to the beacon system. To ameliorate this, the sensor head also has affixed 
to it a digital compass [22] that at every sensor location gives the three Euler angles, in the yaw- 
pitch-roll convention [23], in a coordinate system defined by magnetic North and the acceleration 
of gravity. 

The connection between the compass and beacon coordinate systems can be determined in a 
least-squares sense from the fact that the direction of the dipole moment (2 in a reference frame 
fixed to the sensor head) is unique. We thus find 


Ryo = [2 22 ... 2M] [eh 22... ayy! (4.1.14) 


where the dagger denotes the pseudoinverse[24]. This rotation matrix can be polished further (i.e., 
made “more orthogonal’’) by iterating[25] 


R,+(R;,) 1 
2 2 
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Figure 4.1.4: Integrated (blue) and averaged (green) flux profiles on a grid above an &-directed dipole; from 
left to right we have the z-, x-, and y-components. The red is the difference between the two. 


with Ro being the matrix found from (4.1.14). Very few iterations (rarely more than three) are 
needed. 


4.1.1 Forward models 


4.1.d.(1) The dipole approximation A popular and useful model for UXO discrimination re- 
places the buried target by one or more triaxial point dipoles, either in the frequency domain [7, 26] 
or in the time domain [9]. The dipole model has been found to fit measured data adequately when 
bodies are small enough that the primary field can be assumed uniform along their extent, when 
they are far enough away from the sensor that higher moments vanish, and when their composi- 
tion is homogeneous enough that there are no conflicting signals from different parts of the ob- 
ject [8, 27]. 
The signal picked up by the pth receiver is when the sensor is centered at rj is 
3(m-R)R/R?—-m . 
Su(r;) = Ἵ " ( : vi iy dsRu (4.1.16) 


where R = r—r’, r is any observation point, and r’ is the location of the dipole. The dipole 
moment m = BHP". An example fit to actual UXO data using this model, as well as the extracted 
polarizability elements, can be found in Fig. 4.1.5. 

We can use the results given by this approximation as a starting guess for the more demanding 
SEA approach, which we discuss next. 
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Figure 4.1.5: (Upper) Quadrature part of the flux through the three receiver coils for an H-83 UXO with 
nose tilted 45° down. (Lower) Extracted polarizability elements (dashed lines: inphase part; solid lines: 
quadrature part); it is clearly seen that the object exhibits cylindrical symmetry. 


4.1.d.(2) The Generalized Standardized Excitation Approach The Generalized Standardized Ex- 
citation Approach (GSEA) [14, 15] is simultaneously a generalization and a simplification of the 
usual SEA[11]. For determining a superior set of equivalent sourcesto produce the fundamen- 
tal modal responses, its starting point is the observation that in the EMI regime the fields out- 
side a metallic target are irrotational and can hence be described by a scalar potential that obeys 
the Laplace equation, whose fundamental solution in a prolate spheroidal system (—1 < ἡ < 1, 
1<& <o,0< @ < 27) centered at the scatterer is the superposition [28] 


“ὍΤΕ om 
V(r) = pHod YD bomnba'(M)Fa'(S )Tom(), where Ton) = 7 
m=0n=m p= 
(4.1.17) 


and the P’"(-) are associated Legendre functions. The interfocal distance for a prolate spheroid 
with semiminor and semimajor axes a and b is d = 2\/b? — a?; the surface of the spheroid is the 
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set of points € = & = 2b/d. We choose the scale factor Hp = 1. It is easier and more efficient to 
perform the decomposition on the normal magnetic field over the spheroid, employing the fact that 


i lo d 
Hg (r) =—€-Vw(r) = The oe , where the metric coefficient he = 5 Ὁ Ὁ Ξ 
(4.1.18) 
The orthogonality of the Legendre and trigonometric functions [29] yields 
1 20 
[ane [6 Tom(O) ne He 
τὸν 5 (n+m)! aP™(E) ὃ (4.1.19) 
2 2n+1(n—m)! Of ᾿ς 


where @ = 2 for p= m=0 and αὶ = 1 otherwise. To compute (4.1.19) we use for each transmitter 
coil the exact expression (4.1.2) for the field of a circular current and evaluate the integrals using 
Gauss-Legendre quadrature and the periodic trapezoid rule[30]. The signal in the th receiver 
when the sensor is at rj; is given by 


co co N 
= YY Peel (τ (4.1.20) 
m=O0n=m /=1 
where 
=e 3 (é'-R)) R)/R? -- ξ' 
Dit) df dspy tty: Tom(" : ἀφ΄ (4.1.21 
ΓΙ -Σ Ρ J Ru(r;) Ἀμμ 0 Ρ (φ) 4πῈ} Xl φ ( ) 


with R; = r -- ΤΊ] and the position vector on the /th ring 


: 7 d 
r) = χι(ᾶὰ cosg/ +9 sing’) +2 €, and χι = SV =nF ee - 1). (4.1.22) 


In other words, instead of the rings of charge considered before [11, 13] we use rings of radially 
directed dipoles. This is a simplification of the original SEA in that it is not necessary to force the 
total charge to vanish in order to comply with the Maxwell equations. Equivalently, 


2π ἔ' χΒ, 
b mn (rj ni) f dl ᾿ T, m - d 4.1.23 
ΓΙ -Σ Pp 1) Ru(r;) Ru 0 Pp (φ' ) 4 “ARR? XI φ΄. ( ) 


Equation (4.1.20) can be turned into a linear system by finding an appropriate invertible map- 
ping (/,mn) ok 


T 
S'=[Si(r) --- S.(rv,) Shri) --- δι ν,)]}. (4.1.24) 
and 
Thidti) Ὁ Ty ltr) τ En, tina 1) 
I7 ro ἧς ro | r2 
Γ-- ial ) nt 0"2) ) . (4.1.25) 
Py (tw) +++ Τα αν) Nrskinax ENp) 
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Once we have this, we minimize the discrepancy to obtain 


min 5(Fp—S’)'(Fp—S’). (4.1.26) 
ρ 


This is a highly ill-conditioned system that requires regularization. We have employed two regu- 
larization methods: the truncated Singular Value Decomposition (TSVD) [31] and the Tikhonov 
procedure [24]. 


4.1.e Results 


We start by stressing the necessity of integrating over the receivers when calculating responses 
for the case of monostatic measurements. We consider a sphere. It is well known[32] that a 
sphere of conductivity σ΄, relative permeability U,, and radius a placed in a uniform primary field 
behaves like a dipole of moment m = BH"; for small spheres we can assume H? from (4.1.2) to 


be approximately uniform over the target size, and the polarizability matrix B = BI, where 


2u, + 1)(kacothka— 1) — (ka)? 
TS TEER ΣΝ annie ee — 


B =2na> 


Equation (4.1.6), which is valid only right above a target, still applies in this case and lets us 
check the correctness of our results. The differences are dramatic when the target is very close to 
the sensor, as Figure 4.1.6 illustrates. In this case the depth is 10 cm, and the sphere is below the 
origin. Equation (4.1.6) is satisfied to 15 figures. It is evident that in the interest of realism we 
must integrate over the receiver areas when modeling responses. 

We now turn to measured data taken over actual spheres at three different sensor heights (16, 
21, and 26 cm above the surface of the sphere). We use two of the data files to extract the dipole 
moments of the rings and use these to predict the field in the unseen instance. We have used 
two methods to regularize (4.1.26): 1) We first perform Tikhonov regularization [11, 24] using 
the L-curve method [31] to determine the regularization parameter. 2) We also computed the 
Singular Value Decomposition (SVD) of F'T, truncated it using the criterion suggested by the 
Picard plot [31], and used the truncated version to solve the normal equation (4.1.26). 


4.1. Conclusion 


In this section we have presented a handheld frequency-domain sensor whose wide range of oper- 
ating frequencies, three-dimensional vector receptivity, and innovative positioning system make it 
a strong candidate for providing robust and reliable discrimination of UXO. 

After introducing the GEM-39~ and seeing the importance of obeying Faraday’s law as ex- 
actly as possible when modeling the signal collected by each receiver, we proceeded to use the 
point dipole model and the Generalized Standardized Excitation Approach to model synthetic and 
measured data. 
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Figure 4.1.6: Integrated (blue) and averaged (green) flux profiles on a grid above a sphere; from left to right 
we have the z-, x-, and y-components. The red is the difference between the two. 
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Figure 4.1.7: (Upper) “Training” field used to extract the GSEA dipole moments of an aluminum sphere; 
(Lower) Picard plot showing the singular values of ΓΊΓ. 
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4.2 Soil Response of the GEM-39* 


Our April 2007 White paper for the GEM-3Y* and the MPV instruments details the theory behind 
soil responses for both instruments. In this section, we will discuss some of the data we have taken 
from local soils with the GEM-39* and the advantages afforded by the transverse receiver coils 
is their relative immunity to contamination from geological noise due to ground response. 

We arranged to test some soils at a local gravel pit, Lebanon Crushed Stone. The main soil 
we tested was known as ’sump sand” and it exhibited a large magnetic response in the frequency 
domain to the GEM-39*. We acquired height profiles at 2 in. increments up to 36 in. over 
this soil. Our measurement set up can be seen in Fig. 4.2.1. We also took high profiles over the 
soil after burying a steel sphere northwest of the measurement position (see Fig. 4.2.2). Height 
profiles at different frequencies after background subtraction are shown in Fig. 4.2.3. As we have 
reported elsewhere, the inphase part of the GEM-39~ data displays a trend of increasing signal 
with decreasing height until a maximum at around 5 to 10 cm above the ground. The quadrature 
component of the heights profile, on the other hand, displays only a weak trend if at all. 


Figure 4.2.1: Setup of GEM-3Y* at Lebanon Crushed Stone over “sump sand”. 


Because the primary field is oriented downward in the —Z direction, the principal response from 
the soil is oriented in 2 with only a secondary response in the transverse receivers (see Fig. 4.2.4). 
This effect can be seen by comparing the data from the height profiles with and without this sphere 
present (see Figs. 4.2.5 and 4.2.6). With whatever variations in magnitude, the isotropic response 
from the sphere target should be the same in all vector components. Notice in Fig. 4.2.6 top, 
that the apparent H, response is shifted downward (c.f. bottom plot) as it also contains the soil 
response. These figures show the strong broadband soil response in the Z receiver coil and the 
almost negligible response in the transverse receiver coils (* and 9). 

Using the scale factors contained in Fig. 4.2.7 and derived from the comparison of analytical 
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Figure 4.2.2: Burying a steel sphere in the “sump sand”. 


models with actual sphere data (see Sec. 4.4), we were also able to estimate the susceptibility, 7, 
of this soil. A plot of the match between model and data is provided in Fig. 4.2.8. The estimated 
susceptibility of this soil is χ = 0.0032, a very realistic value based on the work of Janet Simms at 
the Standardized Test Sites [33-35]. 

Using this technique, we hope to be able to account for local geological noise at remediation 
sites by acquiring height profile data over local soil which is not contaminated by metallic targets; 
and also by inferring soil effects in the H, component by analyzing the relatively uncontaminated 
horizontal components. 
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Figure 4.2.3: Inphase and quadrature components of GEM-3Y* taken over “sump sand” with height de- 


creasing left to right. 
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Figure 4.2.4: EMI sensor above ἃ halfspace, including a bump on the soil-air interface. 


Figure 4.2.5: GEM-39* setup comparing “sump sand” data to “sump sand”+sphere data. 
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Figure 4.2.6: GEM-3Y* inphase response showing the substantial soil response contaminating the Ης re- 
ceiver but not the ἢν receiver. 
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Alum | Alum | Brass | Brass | Steel Steel Overall | Max. 
16cm | 21cm | 16cm | 21cm | 17.5cm | 22.5cm | Avg. Dev. 
Avg. |Avg. |Avg. |Avg. | Avg. Avg. per from 
Loop | Avg. 
Hx | 3252 |3136 |3174 |3115 |3203 |3067 |3159 | 2.943% 
Hy / 3338 |3195 |3283 [3153 |3312 |3187 |3245 | 2.866% 


Hz [2129 [2120 |2183 [2084 [2160 |2132 |2135 | 2.248% 


Figure 4.2.7: Soil model match to GEM-3Y* data using scale factor. 


Measurement of Soil, 4 inches above Ground 


Figure 4.2.8: Soil model match to GEM-39* data using scale factor. 
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See additional material relative to soil response — addendum to “White Paper’ for MM1537 [36] 
that illustrates the use of the conversion factors for inference of geophysical conclusions, and for 
characterization of instruments relative to UXO discrimination. 
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The problem with the view encouraged by Fig. 4.2.8 is that it does not explain the fact that 
significant soil responses are seen by time domain instruments. To understand the issue, consider 
a FD expression for μι, in terms of magnetic susceptibility 7 as 


M,(@) =1+%(@)=1+4'(@)+ix" (ὦ) (4.2.1) 


where ἱ is the square root of negative one and the angular frequency ὦ is 27 times frequency in 
Hz. The implication of spectral patterns like that in Fig. 4.2.8 is that 


χ' (@) = constant, x" (@) ~ 0 (4.2.2) 


Such a constant, real valued y should produce no response in time domain measurements. Given 
negligible electrical conductivity in a soil, its response would simply be proportional tox’times the 
magnitude of the transmitted field at any given point in time, with no delay. All response would be 
instantaneous. When the transmitted field from a TD instrument is shut off during the recording of 
data, the soil response would also shut off, producing no data other than noise. However, many time 
domain measurements do indeed show some soil response in more or less orderly decay patterns. 
Evidently some kind of relaxation effect takes place, beyond the instantaneous response. 

Overall, one might hypothesize a soil response in the form of the schematic in Fig. 4.2.9, given 
a downward step input in the primary field at t = 0. 


Schematic of ground response 
sequence in time, with Tx shutoff at t = 0 


u (t) 


Normalized Time 


Figure 4.2.9: Hypothesized relaxation type time response for real soil permeabilities. 


The instantaneous response of the system is from its initial value of u, to u;. This response per- 
sists for all t>0. It corresponds to the flat real-valued (i.e. inphase) impulse response in Fig. 4.2.9. 
On top of that is a relaxation to the ultimate response to Ug. 

The only way to explain measurable TD instrument response to the same soil as in Fig. 4.2.8 is 
to hypothesize that a relaxation type spectrum (see MM1537 white paper) lies within the spectrum 
shown in the Fig. 4.2.8, but that it is overshadowed by the larger instantaneous response, i.e. it is 
below the resolution of this data. 
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To investigate this, we also computed conversion factors for the TD MPV device (see Sec. 4.4). 
If the conversion/calibration factors for both the GEM and the MPV are correct, then when the 
primary fields are scaled so that a sphere is illuminated by the same intensity from each device, 
each instrument will report the same magnetic response at the same distance. Figure 4.2.10 below 
shows that this is the case, when the GEM-3 spectrum is converted to a TD response via Wait’s 
analytical solution [37, 38]. 


1011 ‘ini Sphere TD response 
Γ === (Converted sphere FD response 


Converted Response 
(A’m of converted response)/(A/m of incident field) 


107 


Figure 4.2.10: Measured GEM-3 (“converted”) and MPV responses from a sphere, observed at the same 
distance, when each response has been scaled additionally by the primary field from each instrument. 


The implication of all this is apparent if one applies the same kind of scaling as in Fig. 4.2.10 
to measured soil responses. For a CRREL backyard soil with FD spectra like that in Fig. 4.2.8 
above, we also measured the MPV response. We then converted the GEM-3 soil response to a 
TD response, retaining the magnitude indicated by our conversion factor, with observation points 
about the same distance from the soil. When the GEM-3 and MPV soil responses are further scaled 
so that the illumination of the soil surface is the same magnitude for both instruments, one obtains 
Fig. 4.2.11 below. 

That is, Fig. 4.2.11 shows approximately the same kind of thing as Fig. 4.2.10. The latter shows 
that, given our conversion factors, we can illustrate that a given metal object (sphere) produces the 
same response to the same to the same intensity of illumination, when observed from the same 
distance by different instruments. The reason that the similarly scaled soil responses in Fig. 4.2.11 
are different is that the TD instrument only sees part of the soil response, namely the relaxation 
part. The instantaneous part has passed by the time the first time gate falls in the MPV. On the 
other hand, the GEM-3, and for that matter any other frequency domain instrument, is dominated 
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10 =—=|Vieasured FD response 
===Measured TD response 


Converted Response 
(A’m of converted response)/(A/m of incident field) 


-0.01 -0.008 -0.006 -0.004 -0.002 τὴ 0.002 0.004 0.006 0.008 0.01 
ime 


Figure 4.2.11: Comparison of measured responses of the soil to the GEM-3 and from the MPV, from about 
the same observation point, when the illumination of the soil surface has been scaled to be the same for each 
instrument. The red GEM-3 response, converted to TD, is flat before the step shut off at time zero, at which 
point it theoretically goes to zero, instantaneously. The blue MPV data start some finite time after shutoff 
(~ 10-4 9) and decay thereafter. 


by the much larger (> 3 OM) instantaneous response, that appears at all frequencies. 
With all this in mind, the most general representation of soil susceptibility may be something 
like: 


X(@) = χοι + χο (a) (4.2.3) 
where χοι and 7,2 are constants and R(@) is something like the commonly used expression: 
1 iWT™+ 1 
R(@) =1- i 4.2.4 
νῷ; [π(το τι) n( Bae) ( ) 


where τὶ and 72 are chosen relaxation time limits. The parametery,; would account for the larger 
flat response in the figs above. 7,2would be the separate magnitude of the relaxation response. In 
this view, all the FD data treated here reflect the fact that χοῦ < Yo. 

A little more explanation in terms of relaxation times. The schematic below in Fig. 4.2.12 
shows hypothetical relaxation responses by different constituents within some soil. 

The commonly used expression in (4.2.4) corresponds essentially to the collection of blue 
curves. That model assumes that there is a continuous distribution of relaxation times, evenly 
distributed in log space, and that all the curves converge to the same magnitude of contribution. 
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Initial value before shutoff 


Slower relaxation 
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("instantaneous") 
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Figure 4.2.12: Hypothetical relaxation curves in (log) time for different constituents of some soil. 


This has been said to model the TD responses effectively for many soils encountered in UXO 
studies [39-41]. What our research and conversion factor applications show is that an additional 
component with very short relaxation time is also typically present, converging to a much larger 
magnitude response (red curve/line). Treating the fast response as effectively instantaneous simply 
means adding the term 7 in (4.2.3). 

It has been suggested that one might cover all the effects mentioned here while retaining the 
older soil response model, i.e. (4.2.3)-(4.2.4) with χοι = 0, by stretching the relaxation time limits, 
in the very least to encompass virtually instantaneous responses. However this is likely to produce 
a muddled formulation, with limited generality. It would imply a continuous distribution of relax- 
ation times between the “instantaneous” one and all the others, with equalization of magnitudes. 
Basically, these “cures” of the model’s shortcomings would constitute unphysical tweakings of 
the parameters or, in the least, an unsupported flight. While they might succeed in matching an 
individual FD spectral curve, the results would be unlikely to generalize successfully. That is, 
application of the resultant 7 to other instruments or circumstances would likely fail. 

CONCLUSIONS/ SIGNIFICANCE: 


1. Application of the scaling/calibration/conversion factors so we can compare FD and TD re- 
sponses in the same terms provides a fundamental insight into how different kinds of in- 
struments respond to significantly different aspects of soil magnetic susceptibility. For the 
various backyard soils examined thus far, results suggest that the slow relaxation TD re- 
sponse is objectively much smaller than the FD response, which is dominated by effectively 
instantaneous response components. 


2. The calibration factors that produced the evaluation of relative magnitude of FD/TD re- 
sponse, as in Fig. 4.2.10, can easily be applied further to calculate the factors in an expanded 
soil magnetic susceptibility model, as per (4.2.3). This information is then potentially trans- 
ferable to other EMI instruments of whatever type. As the material above shows, existing 
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models are not transferable. Heretofore, FD data implied no TD response; TD data implied a 
much weaker, slower response than was seen in the FD. The expanded model covers all these 
bases. 


3. Item #1 suggests that FD instruments will be at a disadvantage for discrimination in instances 
where significant ground responses cannot readily be subtracted as background. The GEM- 
3D+ counteracts this disadvantage by recording horizontal field components that are largely 
immune to ground response and that can help distinguish the contribution of target vs ground 
in typically stronger vertical components. 


4.3 Estimating magnetic susceptibility from EMI data 
4.3.a Executive Summary 


Studies have showed that magnetically susceptible soils significantly affect EMI sensor perfor- 
mance, which in return reduce the sensors discrimination capabilities. In order to improve EMI 
sensors detection and discrimination performances first soil’s magnetic susceptibility needs to be 
estimated, and then the soils EMI responses have to be taken into account during geophysical data 
inversion procedure. Until now the soil’s magnetic susceptibility is determined using a tiny amount 
(up to 15 mg) of soil’s probe. This approach in many cases does not represent effective magnetic 
susceptibility that affects on the EMI sensors performances. This paper presents an approach for 
estimating soil’s magnetic susceptibility from low frequency electromagnetic induction data and it 
is designed namely for the GeoPhex frequency domain GEM-3 sensor. In addition, a numerical 
code called the method auxiliary sources (MAS) is employed for establishing relation between 
magnetically susceptible soil’s surface statistics and EMI scattered field. Using the MAS code 
EMI scatterings are studied for magnetically susceptible soils with two types of surfaces: body of 
revolution (BOR) and 3D rough surface. To demonstrate applicability of the technique first the 
magnetic susceptibility is inverted from frequency domain data that were collected at Cold Re- 
gions Research and Engineering Laboratory’s test-stand site. Then, several numerical results are 
presented to demonstrate the relation between surface roughness statistic and EMI scattered fields. 


4.3.b Introduction 


Active EMI sensors induce magnetization via their primary fields. The static magnetic sensors 
are sensitive to the presence of magnetic minerals because the earth’s magnetic field creates an 
induced magnetization in them. This anomalous magnetic field produced by the soil magneti- 
zation will be superimposed on the magnetic anomalies produced by UXO’s and other metallic 
objects and will therefore introduce adverse effects into the discrimination problem. There are 
three magnetic effects that impact the magnetic and electromagnetic characteristics of the subsur- 
face: (4.3.1) induced magnetization, (4.3.2) viscous remanent magnetization (VRM); and (4.3.3) 
permanent/remanent magnetization 

The magnetically susceptible and heterogeneous soils are a major problem for the near field 
EMI sensing technologies. Studies [42—51] showed that, magnetically susceptible soils can pro- 
duce electromagnetic anomalies of the same magnitude as buried metallic targets. Even that in 
those studies, it is assumed that the spatial distribution of magnetic anomalies are constant, similar 
to a half space, and that this response can be subtracted from measured data, the practice shows that 
the magnetic soils significantly degrades performances and applicability of the EMI technologies 
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for UXO discrimination. Since, for discriminating between UXO and non-UXO targets accurately 
a dense measurement grid is required. Under these conditions, variations in sensor height and ori- 
entation, as well as small variations in the surface topography, can produce anomalies similar to 
those from UXO [44]. However, in most existing approaches to buried metallic object discrimina- 
tion, the object of interest is assumed to be embedded in a medium with the same electromagnetic 
properties as free-space. Any influence of the background medium assumed to have been removed 
by pre-processing or filtering before the data are submitted to an inversion routine. Until now, 
techniques for calculating the background effects, which include the half layered-space earth and 
an equivalent dipole layer models that attempts to model the background geology, have almost 
exclusively neglected soils surface roughness and these techniques are limited to the use of one 
dimensional or layered Earth models. Such models are incapable of capturing the true process 
of geomagnetic induction that is associated with the surface roughness, as it is the case in real 
field environment. Therefore, one must characterize soil response so that it can be predicted, or 
so that it can at least be accounted for during data processing. The characterizing of soil’s EMI 
response requires knowledge of its magnetic susceptibility. To date the soil’s magnetic suscepti- 
bility is measured using a very small amount (up to 15 mg) of soil’s probe, which in many cases 
does not represent effective magnetic susceptibility that affects on the state of the art EMI sensors 
performances. 

This section investigates EMI scattering from magnetically susceptible random rough surfaces 
for establishing statistical relations between surface and signal characteristics. The study is done 
for the handheld GEM-3 sensor. First, to understand interaction between magnetically susceptible 
soil and the sensor here, an EMI magneto quasi-static (MQS) problem is solved is solved using 
a numerical approach called the method of auxiliary sources (MAS) [52-55]. In MQS magnetic 
fields are irrotational, and can thus be represented efficiently using a magnetic scalar potential. 
In the MAS, boundary value problems are solved numerically by representing the electromagnetic 
fields in each domain of the structure under investigation by a finite linear combination of analytical 
solutions of the relevant field equations, corresponding to sources situated at some distance away 
from the boundaries of each domain. The ’auxiliary sources” producing these analytical solutions 
are chosen to be elementary dipoles/charges located on fictitious auxiliary surface(s), usually con- 
forming to but offset slightly from the actual surface(s) of the structure. Enforcement of standard 
electromagnetic boundary conditions at an array of points over the object’s actual surface allows 
us to solve for the auxiliary sources, from which we can immediately express all EM fields in the 
problem. 

The section is organized as follows: In section II, overviews very briefly EMI problem and 
the MAS Galerkin technique for MQS regime, Section III describes an approach for estimating 
soil’s magnetic susceptibility from EMI data, section IV demonstrates several numerical results, 
and finally section V, presents conclusions and discusions. 


4.3.c EMI Scattering for Magentic Rough Soil 


To illustrate the MAS for MQS problem, let us consider a permeable 3D magnetically susceptible 
rough surface. The surface is illuminated with an electromagnetic field produced by a GEM- 
3 sensor (Fig. 4.3.1). The computational region is divided into two regions: the region above 
the surface is referred to as Region 0 with free-space electromagnetic parameters, and Region 
1 corresponding to the volume under the rough surface with €;=1 but with a nonzero magnetic 
susceptibility. As is well established in the MQS regime, the displacement currents are negligible. 
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The primary field induces magnetic dipoles within the Region 1, which in turn produce secondary 
(scattered) fields. The magnetic field inside each region can be expressed as 


Hq(?) = —VPa(7) (4.3.1) 
where 
Ba (F) ore ! ἢ =A (4.3.2) 
χρο 


Here Ψα() is the fundamental solution of the Poison equation, @=0,1, μα is the permeability 


GEM-3 


Figure 4.3.1: 3-D rough surface used for computational model. 


of region @, and ρα is the πῇ! unknown magnetic charge source coefficient distributed on the a 


auxiliary surfaces. These sources reside on an auxiliary surface slightly outside and approximately 
conforming to the real surface. We compute the total field inside each ἃ = 0,1 region by assuming 
that these sources radiate into an infinite homogeneous medium with the properties of the corre- 
sponding region. 

The electromagnetic fields at each interface between two regions are connected by boundary 
conditions. Specifically, the magnetic scalar potential and the normal component of the magnetic 
flux density must be continuous across the real surface: 


[ἢ x (Wier — wer )] = 0 (4.3.3) 


᾿Ξ (By πρίοἵ 


2B) = ὃ (4.3.4) 


where fi is the unit normal vector on the interface between the two regions. For achieving high 
accuracy and saving computational resources, the boundary value problem is solved using the MAS 
together with the standard Galerkin technique. To demonstration the MAS-Galerkin approach, let 
us rewrite the boundary value problem in the following compact form 


Lp=vy (4.3.5) 
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Here L is an integral-differential operator, p is auxiliary magnetic charge density and y is the right 
hand side of (4.3.2). Let a weighing function w and (f,w) denote the inner product of two scalar 
functions fand w, usually defined as 


(fw) = / ἐμὴν (4.3.6) 


here ds’ is a surface element. Then, by using Galerkin’s method, one can rewrite equation (4.3.6) 
as 
(Lp,w) = (W,w) (4.3.7) 


upon discretization of equation (4.3.7) a linear system of equations for unknown p coefficients 
results. Once the amplitudes of the p coefficients are determined, the field at all points in both 
regions can be readily computed. 


Figure 4.3.2: Experimental setup of GEM-3 sensor at various heights above the ground. 


4.3.d Estimation of Soil’s Magnetic Susceptibility from GEM-3 Measurement Data 


The GEM-3 sensor consists of two transmitter loops. Based on the standard expression for the field 
from a wire loop in the magneto-quasistatic regime, the z component of the magnetic field radiated 
by the GEM-3 sensor along its z rotational axis can be expressed as: 
1 2 2 2 
I; ay + 1; a5 
3/2 3/2 
2 (a? +22) Pp 2 (a3 +27) f 


(4.3.8) 


Ζ --- 


According to image theory, for ἃ GEM-3 sensor placed at an elevation z above ἃ permeable half 
space, the magnetic field produced at z by the ground can be written as 


η, -- ἵει -ὶ Ἰραῖ : = : (4.3.9) 
1+ [gr 2 (a2 :. 4.2) (a? + 1622) *” 
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where 11 = I,,andI? = —I,/2 has been assumed. So, if we know the J, current and the a, radius 
then we can calculate the permeability of ground [,,. However, different EMI sensors record 
data in different units. For example, the GEM sensor produces field in “ppm” units relative to a 
separate reference receiver coil. Given the known geometry and assumed relation of the currents 
in the GEM, the actual scattering field from any target (metal or ground) with o and yp, can be 
written in the following form 


H = IF (0,L,geometry) = Haata/P (4.3.10) 


Tym 15 the unknown current and p is some unknown scaling factor (accounting for amplification, 
filtering, the size of the reference coil, etc.) that produces the GEM’s readings in terms of “ppm”. 
Once the scaling factor is known, then the soil’s magnetic susceptibility can be estimated as 


X = Uer-l = — (4.3.11) 


here 


Haata 1ραῇ 1 1 
ys ae = (4.3.12) 
PB 2 (aj -Ε 422) ie (aj + 16:2) = 


Using methods explained in Sec. 4.4 the scaling factor p was determined in [56] for the GEM- 
3 sensor with J, = 10A current to be p = 4.32- 10°. This value was determined from detailed 
measured data for a sphere. 
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Figure 4.3.3: Mean value of the predicted magnetic susceptibility versus frequency. 
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Figure 4.3.4: 1-D Rough surface. In meters. 
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Figure 4.3.5: Scattered magnetic field for 1-D rough surface. 
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4.3.e Results 


Experimental data were collected using GEM-3 sensor at CRREL facilities. To estimate soil’s 
magnetic susceptibility the measurement were taken over a particular spot as a function of sensor 
elevation. For this the “blue rig” (see Fig. 4.3.2 left) was used. For the test 18 shots were taken as 
follows: three on the top of the rig, as in the picture above left; 9 descending from the top to within 
10 cm of the base; then the rig was removed and 3 shots were taken at successive holes about 10 
cm apart on a board placed on the ground where the rig had been (Fig. 4.3.2 right); then the blue 
rig was put back in position and 3 shots were again taken on the top. The mean value (respect to 
the evaluation) of predicted soil 

susceptibility real and imaginary parts for CRREL soil versus frequency is shown on Fig. 4.3.3. 
The results show that the soil has a complex magnetic susceptibility. The amplitude of the real part 
is around 0.38 x 10~, and it is decreasing as frequency increases, as expected. The imaginary part 
also shows very small but frequency dependent trends, indicating that the present soil will affect on 
the time domain EMI sensors performances as well. We intend to deploy the time domain EM-63 
sensor on the same soil in near future, and outcome investigations will be reported in subsequent 
publications. 

Once the soil’s magnetic susceptibility was determined then the MAS code was used for GEM- 
3 sensor for establishing statistical relation between surface roughness and EMI scattered field 
for magnetic soil. First BOR rough surfaces were generated using 1d code [57] for generating 
14 BOR rough surface. A snap shot of the surface in a cylindrical coordinate system is depicted 
on Fig. 4.3.4. On this figure vertical and horizontal axis coincide with the z and p axis in the 
cylindrical coordinate system, respectively. 1D line on Fig. 4.3.4 was rotated around the vertical 
axis for generating the BOR surface. The GEM-3 sensor that is the excitation source was placed 
at 15 cm distance from the origin (0,0 point on Fig. 4.3.4). 

For each fixed correlation length and RMS heights 100 realizations were simulated. One of 
such runs shown on Fig. 4.3.5. The standard deviation and mean are shown on Fig. 4.3.6. These 
results show that as correlation length and RMS height increases the EMI response from magnetic 
soil increases as expected. The results also show that for high correlation length (surface is becom- 
ing more flat) and for high RMS height (sensor is close to the surface), the standard deviation of 
the EMI response decreases slightly. 

Finally, we studied EMI scattering for 3D rough surfaces. First, we construct 3d rough surfaces. 
One of rough surfaces is shown on Fig. 4.3.6. The GEM-3 sensor was placed 25 cm above the 
surface. Using the MAS-Galerkin technique the full EMI scattering problem was solved. For each 
fixed correlation length and RMC height EMI scatted filed was computed for 100 realizations. 
Figure 4.3.8 shows the scattered magnetic field versus realization number. The standard deviation 
and mean values are depicted on Fig. 4.3.6 for 3d rough surface, with fixed correlation length. The 
results show that both standard deviation and mean values increase as RMS height of the 3d rough 
surfaces increases. 


4.3. Conclusion 


In this section, the EMI secondary response from a magnetically susceptible soil with rough sur- 
faces were investigated and analyzed. The combined MAS —Galerkin technique is employed for 
investigating EMI scattering from 3d-rough surfaces. The relation between surface statistics and 
EMI scattered field is investigated and presented. In future work the relation between sensors el- 
evation, orientation and surface roughness statistic will be studied and reported in a subsequent 
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Figure 4.3.6: Standard deviation and mean value for the BOR rough surfaces versus RMS height and corre- 
lation length. 
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Figure 4.3.7: 3-D random rough surface. 
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Figure 4.3.8: Scattered EMI field from 3d-rough surface. 
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Figure 4.3.9: Standard deviation and mean value for the 3d rough surfaces versus RMS height at fixed 
surface correlation length. 
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4.4 Absolute Calibration of EMI Measurements and Application to Soil Magnetic Suscep- 
tibility Inference 


This section presents a method for determining a conversion factor for electromagnetic induction 
(EMI) survey measurements so that fields received in an instrument’s native units can be translated 
into standard (“absolute”) EM units. The method is non-invasive, requiring only a rudimentary 
characterization of the instrument, such as loop geometry and the basic nature of the output such 
as the magnetic field, differencing of magnetic field, or time derivative of magnetic field. It does 
not require knowledge of the device’s internal electrical engineering features, such as amplifica- 
tion, filtering, or transfer functions between components. The technique proceeds by comparing 
data from controlled measurements to model results and thus allows one to do similar comparisons 
henceforth. As an example application and reasonableness check, we use a particular frequency 
domain (FD) instrument to infer soil magnetic susceptibility value in-situ directly from survey data. 
The same methodology is applied to a new time domain (TD) instrument. This allows validation 
of the basic scaling methodology via benchmark cases and also illustrates its transferability. Con- 
version of data from both instrument types into standard units, for measurements from the same 
soil, illuminates the contrasting nature of the soil responses to TD and FD sensors. It also points 
to the fundamentally different magnitudes of different components of soil magnetic susceptibility 
in the soil studied, relating to instantaneous response and relaxation response. 

We acquired data over canonical targets and over some different soil types when we first ob- 
tained the GEM-39" . Using these data we were able to extract scaling factors for the three receiver 
coils on the GEM-3Y~. We were also able to estimate the susceptibility of the soils we measured 
directly from the measured data. 

Data Collected 


e Data taken over Spheres and other canonical targets 


e Data collected over 2 different soil types 


4.4.a Introduction and Background 


To understand the needed for an ability to convert a sensor’s output into standardized units, it 
is helpful to view the UXO problem from an inversion perspective. Inversion studies use the 
secondary field to infer the target’s properties. Possibly the most basic example is the estimation 
of soil properties, particularly broadband magnetic susceptibility. EMI responses caused by this 
property of soil can be a major problem in surveys seeking to discriminate UXO [45, 48, 50, 58, 
59]. The ability to estimate susceptibility, in situ, directly from sensor data would provide valuable 
soil characteristic data in an extremely efficient manner. A commonly used method of inversion 
is to match the measured secondary response with an analytical model’s response of homogenous 
or layered half-space media. The susceptibility of the soil will correspond to the values of the 
model’s susceptibility parameters that produce the best match. In many prior studies of this nature, 
the secondary fields are reported in units of parts per million (ppm) of the primary field intensity 
[60-65]. 

Examining the past work reveals two key aspects which differentiate our work. First, all in- 
struments in these prior studies use receiver and transmitter coils consisting of laterally separated 
horizontal current loops. To our knowledge, no similar study has been made with more compact, 
complex instruments where the transmitting and receiving coils may be coaxial. Secondly, more 
significantly, the prior studies assumed that the results reported by their instruments were truly 
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ppm relative to a specific primary field sample. However, modern EMI detectors include filters, 
amplifiers, etc, attached to particular components, the details of which a general user has no way 
of knowing. Therefore the user has no real way of knowing the relation between actual transmitted 
and received fields in A/m or Teslas and the actual units of measurements. 

This situation argues for a simple, fast, non-invasive method of both assessing the performance 
and calibrating the measurements of instruments that do not report data in known, universal units 
to produce a match with the secondary fields predicted by a model. The latter is required to 
produce scientific inferences, such as estimates of ground magnetic susceptibility using the method 
described earlier. Thus, the calibration essentially amounts to translating reported ppm data to be 
consistent with the inversion researcher’s model. 

The need for this translation factor within UXO inversion studies can be seen in the research 
to discriminate UXO from clutter. Machine learning techniques can make use of the modeled 
responses of a wide range of metallic objects to train a systematic discriminator that separates 
UXO from clutter objects on the basis of their EMI responses [66, 67]. To classify objects with 
measured data, previous studies frequently normalize all measured and simulated data [67] to 
ensure consistency. This normalization entails a loss of valuable information about the magnitude 
of the response. For example, a very small piece of clutter may likely produce a weaker signal 
than a large UXO if both are equidistant from the sensor. Normalization would erase such basic 
distinction. However, if one were able to convert data into units consistent with the standardized 
units of the modeled data—in other words, to ensure that the modeled responses and measured 
responses are directly comparable—no information would be lost. Therefore, the recovery of this 
conversion relationship is a prerequisite for many types of advanced inversion studies. 


4.4.b Test Instrumentation 


As a handy example of an ultra-wideband frequency domain (FD) EMI instrument, here we use 
the GEM-3 sensor, manufactured by Geophex [4, 64] and used in inversion research [68, 69]. Our 
choice of the GEM-3 is largely due to availability and relative ease of use. As mentioned earlier, 
this section illustrates a general method of calibration that may be applied to many EMI instru- 
ments. This analysis of the GEM-3 instrument can be directly applied to any non-adaptive EMI 
instrument which reports measurements proportional the transmitted primary field or which uti- 
lizes a fixed but possibly unknown transmitted field strength. Following a detailed treatment of the 
method in application to the GEM-3, we then apply it to a new time domain (TD) instrument, the 
Man-portable Vector Time Domain Electromagnetic Sensor (MPV). This device has just recently 
been developed in an ongoing research project [70]. 

To determine, through non-invasive means, a conversion factor that can translate our particular 
GEM-3 instrument’s record of the secondary field into standard units, we compare the actual re- 
ceived fields from a canonical object to the fields determined by a model. Thereby, this analysis 
permits one to acquire the ability to evaluate the sensor’s data in standardized quantities in all other 
future applications. Furthermore, we examine the significance of the discrete size of the GEM-3 
receiver loop in the isolation of the conversion factor. Ultimately, with proper calibration, data in 
the GEM-3’s native units from magnetically responsive soil is used directly to estimate broadband 
magnetic susceptibility of the soil. 

As an additional, comparative test we utilized a new time domain (TD) instrument designed and 
manufactured under the direction of one of the authors [70]. This time domain instrument consists 
of two coaxial coils of 37.5cm radius that generate a step down function. From 46 microseconds 
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Figure 4.4.1: Flow chart of instrument’s process, to be paralleled by the modeling sequence. If normalization 
is performed by the instrument then the lower, shaded set of boxes may be deleted in the model and an 
arbitrary (e.g., unit) current assumed. 


to about 20 milliseconds after shut off, the step down response is then measured at a receiver loop 
located at the center of the transmitter coils. 


4.4.c Recovery of the Conversion Factor 


A given transmitter (Tx) current in any instrument will produce a particular magnetic field at a 
target’s location (e.g. in A/m), which in turn will produce a particular secondary field value at the 
instrument’s receiver (Rx). In the case of the GEM-3, the Rx signal is normalized, frequency by 
frequency, by a sample of the primary field in a secondary Rx coil. Thus, in effect, the instrument 
provides a signal at each frequency that is proportional to received field per unit Tx current, so we 
need not know the specific Tx current to model what occurs. We do require, however, knowledge 
of the geometry of the transmitter and receiver loops. 

Through our analysis, we can confirm if the proportionality between the actual received field 
and the sensor measurement signal is truly fixed. While the signal may be selectively amplified 
or filtered, this manipulation is always done in the same way, without adaptive measures which 
would alter that relation. Therefore, the recovered conversion factor should have no dependence 
on measurement position or object type. 

We can readily calculate the response of standard objects to fields produced by arbitrary as- 
sumed fixed Tx currents in the known Rx geometry, and thus determine the proportionality be- 
tween measured and actual fields. The standard objects consist here of three spheres made of 
aluminum, brass, and steel. The underlying physical process by which the instrument produces 
signals is presented symbolically in the flow chart of Fig. 4.4.1 and equations that follow. 

Production of the instrument’s signal, in the case of the GEM-3, is 

Ix Fry: Ayys- S -As_.Rx Ὁ Fry 


G Ξ 4.4.1 
(ppm) Tyy: Fry: ΠΤ ΚΕΡῚ ἔΚΕΕ ( ) 


where 1γχ is the current in the transmitter loops. Fr, is what the instrument does with this current 
in producing the Tx field. In the case of the GEM-3 this contains little besides the geometry of 
the coils, as we may simply assume that the current is continuously on during the measurement 


Benjamin Barrowes, Kevin O’ Neill -68- 2 fT] & he) δ: 
E/) AT DARTMOUTH 


MM-1537, GEM-39*- Final Report 4 MATERIALS AND METHODS: GEM-39Y* SOFTWARE 
4.4 Absolute Calibration of EMI Measurements 


portion of the cycle. Hy, 5 1s a transfer function yielding the field produced by the transmitter at 
the target location. S is the inherent target response to a unit input, while Hs_~,, is the function 
that converts that response to the field at the Rx location. FR, is the process that is performed on 
that field by the receiver, including any associated amplification, filtering, etc. The quantities in 
the denominator are explicable in similar terms, except that REF refers to the reference coil, in a 
fixed position relative to to the Tx coils, where the Tx field is sampled. We do not need to know 
any of its particulars. 
The process producing the actual (modeled) field M at the Rx location can be represented as 


M(A/m) = Im: Fry: Ary ss: S + Hssrx: Fre (4.4.2) 


where J, is the current assumed in the transmitter. Fg, may include integration of the secondary 
field over the area of the Rx coil but, unlike Fr, for the actual instrument, need not contain any other 
functions performed by the sensor’s electronics or processors. Note that both G and M contain the 
same product Fr, Hy, —.5: S -Hs5_,r,, which we consider known (calculable). The ratio of recorded 
data to actual (calculated) field values is thus 


G ( Frx ) I (4.4.3) 


M Fry: ΠῚ rer * PREF In FR. 


The product J, - Fe. is known; we may also normalize by taking J, to be unity and adjusting M 
accordingly. The quantity in parenthesis is a function of the internals of the instrument and is what 
we infer as a lumped factor from measurements of the ratio on the left. When this is done, the 
conversion factor has been obtained. 

This method to recover the conversion factor can also be used for any comparable instrument 
by accounting for its particular loop geometry. Note that, to be fully comparable, an instrument 
must have the signal normalized by Tx field or Tx current or use a known, fixed Tx field or current. 
Also, the GEM-3 data is processed within the instrument so that it is proportional to the Rx field 
value, not to the time derivative thereof. Simple alterations of the method can be brought to bear 
to account for differences in any of these particulars. For example, in our TD MPV instrument, the 
current is shut off prior to measurement; the data is not inherently normalized by the Tx current or 
field; and the output is proportional to the time derivative of B and not to B. This simply means 
that we must heed these facts when comparing the calculated and measured responses from our 
standard objects, making sure that the processes parallel one another by adjusting, respectively, the 
content of Fy,, the denominator in equation (4.4.3) , and the contents of Fr, ἀπά μι. 


4.4.c.(1) Modeling of Metallic Spheres For the frequency domain, we model metallic spheres as 
a special case of spheroids in the magnetoquasistic regime. This analytical solution is well vali- 
dated and produces EMI secondary field predictions for spheroidal objects with any conductivity 
and permeability values and in response to arbitrary excitation [71-73]. The solution provides 
secondary field predictions in terms of A/m. The values can be obtained at any observations point, 
e.g. at the nominal location of an Rx coil or over arrays of quadrature points within its area. The 
three specific spheres that were modeled were all 3 inches in diameter. In addition, we imple- 
mented the time domain step down response of a metallic sphere to a uniform field as given by 
Wait and Spies, 1969[38]. Given the relatively small size of the sphere and the distance between 
the sphere and MPV transmitter loops in all measurements, the incident field is approximately 
uniform. Therefore, Wait and Spies’s solution is a valid model. 
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4.4.c.(2) Modeling of the Primary field The excitation field is generated by an analytical model 
that represents the GEM-3 instrument as a set of idealized wire loops. The transmitting sensor 
head consists of two concentric current loops in series with radii of 20 cm and 11.074 cm. The 
smaller loop is placed inside the larger loop and serves as a bucking coil, approximately nulling the 
primary field at the center of the sensor head, where the receiver coil resides. There are 8 windings 
in the outer and 4 windings in the inner transmitter coil. We assume in the calculations that the 
wires in the transmitter coils carry 10 Amps, which is the case at most frequencies. However, as 
mentioned earlier, this varies somewhat by frequency but because of the normalization by Tx field, 
the results are not affected. 

We use the complete elliptical integral functions of the first and second kind to characterize the 
magnetic vector potential and vector field from a circular current loop in a cylindrical coordinate 
system. This calculation method is well established [19, 74, 75]. For a current loop of radius a and 
normal to the z axis in the z= z, plane and with 10 current flowing in the @girection, the transverse 
and axial fields are, 


= Holo (ς -- ζ0) . r+at (2—Z0)* 
= Holo ~P=e Εἰ 2 Ξε)" 
B62) = sara toe [KYO Gage] (45) 
where . 


4ar 
fe = a a) 


and μο is the permeability of the surrounding medium. The standard cylindrical coordinates of 
radius, rotational angle around the z axis, and elevation are expressed by (7; @, z). K and E are the 
complete elliptical integral functions of the first and second kind, respectively. These equations will 
produce the magnetic field from a current loop at any observation position. The total primary field 
is created from the superposition of two current loops that correspond to the two actual transmitting 
loops in the sensor head. Our implementation of the model agrees with the alternative method of 
using numerical integration and application of Biot-Savart’s Law but has the added benefit of being 
more computationally efficient. 

For the time domain MPV instrument, the analytical TD solution for a sphere’s response to 
a uniform step down function only requires calculation of the equivalent of a static uniform field 
generated by a current loop and observed at the location of the sphere. This can be done through 
the current loop equations provided above where the radius of the transmitter loop now 37.5cm 
and the number of windings is 100 to correspond to the MPV. 


4.4.c.(3) Modeling of the Receiver Loop Given the finite size of the receiver loop, a circular loop 
with a radius of 6cm for the GEM-3, the secondary field will not necessarily be uniform across 
the area of this loop. This is especially true if the object is close to the instrument. Therefore 
all modeled secondary fields are numerically integrated over the loop area to parallel the actual 
secondary field measurement. This integration is calculated by partitioning the loop area into 
concentric rings. Then each ring is again segmented along its length such that the entire area is 
partitioned into curved trapezoids. In our experience, over 300 divisions are needed for consistent 
results. 
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Figure 4.4.2: Block diagram illustrating the search algorithm. 


4.4.c.(4) Measurements Measurements were taken above each metal sphere by placing the GEM- 
3 instrument at points in a vertical grid. This was accomplished by measuring along horizontal 
lines above one another at different elevations. The horizontal point spacing was in 2 cm incre- 
ments from -30 cm to +30 cm from the position above the sphere. Vertical spacing of the lines was 
in 5cm increments from 13 cm to 28 cm above the surface of the sphere. The frequency range was 
from 30 Hz to 47 kHz. This wide range of measurement points both in space and frequency was 
taken to ensure that the recovered scaling factor is consistent: there must not be any variation due 
to object distance, frequency, or choice of target object. In other words, all the data must point to 
a scaling factor of the same value. 


4.4.c.(5) Matching Algorithm We performed a Fibonacci search [76] within the range of possible 
values for the scaling factor and the sphere’s material properties to find the best match between the 
model’s prediction and the data. A Fibonacci search first partitions the search range into two 
halves. The median value within each half is given to the model as input to produce two different 
outputs. Whichever half produces the more accurate output will become the new search space and 
the process repeats. The overview of this algorithm is illustrated in Fig. 4.4.2. Accuracy is gauged 
by utilizing the normalized mean square error as an objective function. This type of simple search 
is relatively fast but usually prone to converging on a solution that is not the best overall match. 
In our investigation that difficulty was not encountered. The matching algorithm is independently 
repeated for the frequency response at each measurement point for each sphere. This independence 
allows us to examine the results for consistency over space and frequency. 

The main difficulty in the search algorithm was that the metallic spheres’ exact permeability 
(μ) and conductivity (o) values were not precisely known. These spheres were purchased off the 
shelf from a third party and not manufactured to any specifications other than that they were of the 
correct size and type of metal. Therefore, to match the model’s EMI prediction with measurements, 
both the scaling factor and uw and o must be simultaneously recovered. As a simplification, the 
aluminum and brass spheres can be assumed to be non-permeable, and thus only their o values 
and the scaling factor need to be recovered. Approximate values of sigma for these materials are 
available from standard textbooks, furnishing a reasonableness check. For the highly permeable 
steel sphere, only the ratio between o and wu affects the secondary response [77]. Therefore, 
relative [1 is assumed to be 100 and the normalized o value is recovered along with the scaling 
factor. Adopting a higher relative value produces no change in the resulting conversion factor. 

Fig. 4.4.3 shows the match between the modeled response and the measurement of the steel 
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Figure 4.4.3: Figure showing the matching between model and measurement of a steel sphere at a single 
grid point. 


sphere. The match does extremely well save for the lowest and highest frequency points. We have 
observed previously that the GEM-3 can have difficulties capturing consistent responses at very 
low and very high EMI frequencies. Therefore, the lowest and highest frequency data are ignored 
during the matching process. Furthermore, we proceeded to recover the scaling factor of the MPV 
instrument in the same fashion and with similarly good results. The matching results for only the 
GEM-3 will be presented to avoid redundancy. 


4.4.d Results 


After performing the Fibonacci search for the material parameters of the spheres and finding the 
best scaling factor for each GEM-3 measurement point, we then averaged the scaling factors for 
all objects recovered at the closest 9 measurement points. We arrived at an overall value of 4.93 x 
10°. Division by this scalar will convert all GEM-3 measurements into the average magnetic field 
that would pass through the receiver coil within our model, given our specific transmitter loop. 
All the data can be converted into a percent error from the averaged conversion factor of 4.93 x 
10°. This variation of the scaling factor can likewise be plotted over space. Shown in Fig. 4.4.4 is 
the variation of the scaling over space for the steel sphere. As the figure demonstrates, the variation 
is minimal with most of the error under 5% error. We see similarly good results for the brass and 
aluminum spheres shown in Figs. 4.4.5 and 4.4.6, respectively. Since the model of the GEM-3 
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Figure 4.4.4: Plot of the variation of the retrieved conversion factor, as a percent error from 4.933 x 10%, for 
the steel sphere as a function of position. 


used a fixed current in the transmitter coils, these results verify that the ppm reported by GEM-3 is 
consistently normalized by some sample of the primary transmitted field. In addition to providing 
a conversion factor, the method also potentially provides a quick way for inversion researchers to 
ascertain if their instruments are functioning correctly with respect to any normalization produced 
internally in the instrument. 

If we can assume the secondary field is relatively constant over the receiver loop, then it would 
be unnecessary to numerically integrate the fields within the loop area. However, comparison 
of the consistency of results with and without such integration shows that significant errors are 
introduced by failure to integrate over the receivers. 


4.4.e Application of the Conversion Factor 


As an example application and reasonableness check, we apply the conversion factor to the prob- 
lem of obtaining an estimate of in situ soil magnetic susceptibility. The solid curves in Fig. 4.4.8 
show the Inphase and Quadrature response of a particular sample of soil queried with the GEM-3 
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Figure 4.4.5: Plot of the variation of the retrieved conversion factor, as a percent error from 4.933 x 10°, for 
the brass sphere as a function of position. 


at 12.5 cm above the surface. The Inphase curve shows virtually no change over frequency, and 
the Quadrature component is negligible. Therefore, we implemented a model of the homogeneous, 
permeable, non-conducting half-space response to the GEM-3 excitation by the image method, as 
per Wait, 1985[78]. The output of this model was converted from model units to GEM-3 ppm units 
through multiplication of the recovered conversion factor. Then a d.c. soil susceptibility value, k = 
70x 10~> SI, was found to provide the best match to data as shown in Fig. 4.4.8. This susceptibility 
value is solidly within the range of reasonable soil values as reported by others [33-35]. 

The example above provides primarily an order of magnitude check, which is often all we may 
require in terms of knowledge of soil susceptibility. It also provides an illustration of FD soil 
response, which we can now compare to TD response. By recovering the conversion factor for 
multiple instruments, one can also make meaningful comparison between data from those instru- 
ments, e.g. for soil response, illuminating the different aspects of the susceptibility that dominate 
the responses of the respective sensor types. To this end, the conversion factor of 6810 for the 
TD MPV instrument was obtained in exactly the same manner as for the GEM-3. Results of this 
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Figure 4.4.6: Plot of the variation of the retrieved conversion factor, as a percent error from 4.933 x 10%, for 
the aluminum sphere as a function of position. 


comparison are shown in Sec. 4.2. 


4.4.6 Conclusion and Discussion 


The formulation presented here shows how, by recourse to simple modeling, one may obtain factors 
to convert EMI instrument output to received field values in standard EM units. As all instruments, 
even of the same type, generally have at least slightly different calibration, the method is best ap- 
plied to each individual sensor to be used. In any case, the formulation makes clear the adjustments 
needed to treat different instrument types. Tests reported here reveal that the finite size of receiver 
loops may have to be taken into account to obtain a precise scaling factor, even when those loops 
are small. 

An appeal of the method is that it does not require any detailed knowledge of the internal 
electronics and processing of the instrument under consideration. As long as basic geometry of the 
transmitter coils is known together with a few other rudimentary instrument characteristics, one 
can obtain actual received field values from data by one-time non-invasive calculation. 
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Figure 4.4.7: Plot of the variation of the retrieved conversion factor, as a percent error from 4.243 x 10°, for 
the aluminum sphere as a function of position. No integration was performed over the receiver loop. 


As an example application and reasonableness test, we used our retrieved scaling factor to 
match the modeled response of permeable, susceptible soil to measurements of soil and have found 
reasonable soil susceptibility values. This not only supports the validity of our scaling factor but 
also points to the possibility of characterizing soil from sensor measurements alone, avoiding the 
need to conduct laboratory analysis on core samples outside of in situ conditions. By this means, 
one sensor’s data can yield geophysical information directly that is applicable to other instruments. 
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Figure 4.4.8: Comparison between GEM-3 measured inphase (I) and quadrature (Q) ground response and 
prediction using half-space soil model with an applied scaling factor. 


4.5 Generalized Standardized Excitations Approach 


The generalized standardized excitation approach (GSEA) is presented to enhance UXO discrim- 
ination under realistic field conditions. The GSEA is a fast, numerical, forward model for rep- 
resenting an object’s EMI responses over the entire frequency band from near DC to 100s of 
kHz, including all physical relevant physical considerations, such as near and far field phenomena, 
arbitrary material and geometrical heterogeneity, and internal target interactions. It has been de- 
veloped and tested in both the frequency and time domains for actual UXOs placed in free space. 
The GSEA, which uses magnetic dipoles instead of magnetic charges as responding sources, is 
capable of taking into account the background medium surrounding an object. Given a modeled 
UWB frequency domain (FD) response, the corresponding time domain (TD) response is easily 
obtained by the inverse Fourier transform. Thus the technique is applicable for any FD or TD 
sensor configuration and can treat complex data sets: novel waveforms, multi-axis, vector, or ten- 
sor magnetic or electromagnetic induction data, or any combination of magnetic and EMI data. 
Host media effects are taken into account via appropriate types of Green’s function and equivalent 
dipole sources. Comparisons between simulations and experimental data illustrate that the GSEA 
is a unified approach for reproducing both TD and FD EMI signals for actual UXOs. The EMI 
response from a soil that has a frequency-dependent magnetic susceptibility is studied. The EMI 
responses in both FD and TD domains are analyzed for the model of an actual UXO that is buried 


Benjamin Barrowes, Kevin O’ Neill -77- - ENGINEERING 
&s ΑΚ...) AT DARTMOUTH 


MM-1537, GEM-3 2): -- Final Report 4 MATERIALS AND METHODS: GEM-39* SOFTWARE 
4.5 Generalized Standardized Excitations Approach 


in a magnetically susceptible half space. 


4.5.a Introduction 


Cleanup of unexploded ordnance (UXO) is greatly complicated by the difficulties in distinguishing 
UXO from non-UXO metallic debris, found at most UXO sites. The problem becomes much more 
complicated when signals are contaminated by noise that originates from magnetically susceptible 
and electrically conductive soils [42—48, 50, 51, 79-81]. Until now, in most existing approaches 
to UXO classification, the object of interest is assumed to be placed in a free space [55, 68, 82— 
90]. Any influence of the host medium is considered to be removed by filtering before data are 
submitted to an inversion algorithm. Recent studies show that “geologically hostile” sites cause 
significant problems for magnetometers and electromagnetic induction (EMI) sensors, in terms 
of both decreased probability of detection and increased probability of false alarm. In regions 
of highly magnetic soil, magnetometry and electromagnetic sensors often detect large anomalies 
that are of geologic, rather than of metallic, origin. For instance, [79] documented the problems 
encountered at the former Naval Training Range on Kaho’olawe Island when using EMI sensors. 
During production surveys at the site, approximately 30% of identified anomalies were from false 
positives due to geology, attributed to the strong magnetic viscosity exhibited by the basaltic soils. 

There is also a need for detection and discrimination of UXO in undersea environments [91]. 
In this case, unlike land surveys where the conductivity (< 10~2[S/m]) of non-permeable soils can 
be neglected, it is impossible to neglect EMI responses due to the conductivity of seawater and the 
ocean bottom in marine surveys. All of this leads to a high level of false alarms, which translates 
into an increased workload because each detected anomaly must be treated as if it were an actual 
UXO. Therefore, innovative discrimination techniques that apply to any field condition and that 
reliably, quickly, and accurately distinguish between hazardous UXO and non-hazardous metallic 
items are required. To address this issue here we present a generalized standardized excitation 
approach (GSEA) that is suitable for complex data sets: novel waveforms, multi-axis, vector, 
tensor magnetic or electromagnetic induction data, or any combination of magnetic and EMI data. 

The SEA for objects placed in free space is described in detail in [55, 88, 90]. This work 
extends the SEA for more general cases by using magnetic dipoles as responding sources instead 
of magnetic charges as in [55, 88, 90]. By using magnetic dipoles with the corresponding dyadic 
Green’s function the GSEA becomes applicable for objects placed in a conducting as well as a 
permeable host medium. 

The GSEA can be briefly outlined as follows. For any given object the amplitudes of the re- 
sponding magnetic dipoles are determined and sorted in the universal library for any number of 
basic spheroidal modes. Then, any primary field is decomposed into a set of basis excitations, 
which are then multiplied by appropriate weights (e.g., spheroidal modal decomposition coeffi- 
cients from which we can calculate the target’s complete response just by superposing responses 
of each basic excitation). The key element in the GSEA is to determine the amplitudes of the re- 
sponding sources, which are characteristic only of the object, and are independent of sensor type, 
object location and orientation, and transmitted waveforms. 

There are two ways to determine the amplitudes on the responding magnetic dipoles: (1) using 
measurement data and (2) solving the full 3-D EMI problem in detail. The most straightforward 
way to determine the amplitudes of the responding magnetic dipole sources is to solve a standard 
inverse problem based on the measured data. Obviously, this process requires very good experi- 
mental conditions and a sufficient number of independent measurements of an object of interest in 
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Figure 4.5.1: Problem geometry and reduced set of sources Bp?” distributed along rings on the spheroid 


surface. 


order to reduce the degree of ill-posedness. The ill-posedness makes the solution inaccurate. Re- 
cently, [88] applied such a data-derived approach to extracting the modal response coefficients for 
each candidate by carefully designing the measurements at different distances and orientations in 
free space. However, the accuracy and reliability of the model parameters determined in this way 
may not always be satisfactory due to unavoidable measurement noise and numerical difficulties 
arising from the inherent ill-conditioning of the problem, although a special treatment was applied 
in [88]. In addition, the model parameters were obtained from measurement data with the given 
sensor. Currently all available EMI sensors have certain limitations in both frequency band and 
time. Therefore, the amplitudes of the responding source that are derived from these measured 
data have limitations, they can’t cover all possible EMI and magnetic data, and they could not be 
used to obtain EMI response in time domain for different wave forms and sensors. 

Originally, in [55] and [90], to determine the amplitudes of responding sources, a numerical 
procedure based on the method auxiliary source (MAS) and hybrid MAS thin skin approximation 
(MAS/TSA) was proposed and it is generalized here. In this procedure, by utilizing a full 3-D EMI 
solver, the modal responding coefficients, or strengths of a reduced source set (RSS) [55, 90] are 
determined by employing a physically complete numerical simulation of the object’s response to 
each fundamental excitation mode. The full MAS model-based approach has an advantage over 
the data-derived-based approach, because it is a well-posed EMI problem and it is not dependent 
on measurements. In this approach the RSS can be obtained very accurately for any excitation 
mode and in an ultra-wide band frequency range. Thus, the technique allows users to calculate 
the EM field in both the frequency and time domain and to control the number of input spheroidal 
modes. Here the TD EMI response for a given UXO is calculated directly from FD RSS sources 
just by using the convolution theorem [92-94]. 
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The section is organized as follows: In Sec. 4.5.b, the generalized standardized excitation ap- 
proach is presented, Sec. 4.5.c describes the frequency-dependent magnetic susceptibility model, 
and Sec. 4.5.d shows several experimental and numerical results, demonstrating the applicability 
of the GSEA in both frequency and time domain, and for frequency-dependent magnetic soil. 


4.5.b Generalized Standardized Excitation Approach (GSEA) 


Recently, the SEA has been developed and applied to UXO discrimination [55, 88, 90]. All those 
studies assume that an object is placed in free space. Here, the SEA is generalized to take into 
account conducting and magnetically susceptible host media effects. To illustrate the GSEA, let us 
assume that an object is placed in a background with magnetic permeability μι and conductivity 
σι. Fig. 4.5.1. The object is illuminated by an arbitrarily oriented, time-varying primary magnetic 
field. We surround the object with a fictitious spheroid, which is introduced only as a computational 
aide in the decomposition of the primary magnetic field into fundamental spheroidal modes. We 
choose spheroids because they can assume the general proportions of elongated objects of interest, 
such as UXO, which are also typically bodies of revolution (BOR). Oblate spheroids can also 
be used for flattened shapes. In general, the fictitious surface could be a smooth closed surface, 
as applicable for a related standardized source set approximation described in [55, 90]. On the 
fictitious spheroid given by & = €, (Fig. 4.5.1), the primary magnetic field can be expressed as: 


A” = Σ Σ y Vil sane (4.5.1) 


m=0n=m p= 


The Dymn are coefficients needed to express the primary field and HH’ pmn 18 the pmn mode of the 
primary magnetic field component when bp, =1. The normal component of the primary magnetic 
field on the fictitious spheroid, H: P’ (1,60, @), can be written as 


25 Σ Σ Σ ome n)P, (Eo) Tym ( Φ), (4.5.2) 


m=0n=m p=0 


Hy’ (7,60, 0) 


where (17,6, @) are the standard spheroidal coordinates; d is the spheroid’s interfocal distance, ΡΠ" 
are associated Legendre functions of the first kind, and Tym (@) represents the azimuthal functions 


given by 
= cos (mo), p=0 
Tpm() = { sin(m@), p= 1. (4.5.3) 
By the orthogonality of the associated Legendre functions, the spheroidal expansion coefficients 
Dpmn Can be derived as 


7 2n+1  (n—m)! 7 a 

Dpmn " ymHodP™ (Eq) (n+m)! fee (7) [Ht (n, fo, P)Tpm() “φάη. (4.5.4) 

where y = 2 form = p =O and y= 1 otherwise. The integration in (4.5.4) is evaluated by numerical 
integration. This completes the decomposition of the primary field H ξ "(n, 60,9). 

After the primary magnetic field is decomposed into the pmn spheroidal modes, then the com- 

plete solution for the target to each A?” pmn field is obtained. Since the object is placed in a conduct- 

ing and magnetically susceptible background, the magnetic field in the entire computational space 
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(Regions 1 and 2, Fig. 4.5.1) satisfies Helmholtz’s wave equation and it can be represented with 
magnetic dipoles as 


Ho(7) = | Ga(F7)-PalF)ds’. (4.5.5) 


Here Gy (σ; 7) is the Dyadic Green’s Function, Pg(7’) are amplitudes of magnetic dipoles, and a 
= 1,2 corresponds to region 1 or 2. Here 


Cen aaa 4.56 
all; = ee al’sr),Ga\rsr ~ 4πμαμοκ (4.5.6) 


ko = \/ ωμαμοίωε, + jou) 


where / is the unit dyad, Κα is the wave-number in αἱ region, σὰ and [gare the conductivity 
and relative permeability of a region respectively, and R = [7 -- ΤΊ is the distance between source 
and observation points. We assume that the relative electric permittivities of both regions are 
1; Go (7;7’) is the fundamental solution for the wave equation of the Hertzian magnetic vector 
potential, whereas the exp { j@r} time convention has been implied and suppressed. 

On the surface of the object, total magnetic fields satisfy the following boundary conditions: 


[ax (HS +H?")] = [Ax HD] (4.5.7) 


ἢ - μι (ΕΥ̓ + HP") = ἢ: μ2}} (4.5.8) 


where fis a unit normal vector on the real surface, H’’is the primary magnetic field, H2 is the total 
magnetic field in region 2, and μα is & region’s relative magnetic permeability. 

To determine amplitudes of the magnetic dipoles we have to solve the entire boundary value 
EMI problem for each mode of the primary field just once and then store the amplitudes of the 
responding sources. To do so, one extends each mode to the physical surface within the enclosing 
spheroid and applies the 3-D MAS-TSA method [55, 90]. Finally, the target’s response for each 
primary magnetic field component Bisnis ik is expressed similar to (4.5.5) as: 


Bae) = bmn 3 δὶ (FTP LA (4.5.9) 


where P} Lk is the amplitude of me ΚΙ" auxiliary magnetic dipole, located at the target’s domain 
[24-25], CormespOneing to the A?” prin 


Using Pie k " for each fundamental mode, the total response at any point outside the scatterer can be 
represented as: 


response, and N is the number of auxiliary magnetic sources. 


Ν 
ΒΡ. Σ y Dpmn 2 Gul (, 7) PLL (4.5.10) 


m=0n=m p= k=1 


Thus, after pre-computation of the Pai. ". coefficients for any given object, the EMI scattering 


problem, for any particular 3-D πο ΤῊΝ involving it, breaks down merely to determine the 
spheroidal modal expansion coefficients Dpmn. 

In equation (4.5.10) a substantial number of responding sources Pt Ἵ, are required to represent 
the scattered magnetic field outside the object (including physical surface). Similar to (4.5.10) we 
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can re-express that field quite accurately in terms of a reduced number of sources from the fictitious 
spheroid. The amplitudes of this reduced number of sources [p?’”"”| for each input pmn spheroidal 
mode can be determined by solving a linear system of equations for the normal component of the 
scattered magnetic field as it is shown in [55, 90], and finally the complete secondary magnetic 
field can be represented as 


Nred — 


-Σ Σ Σ bom Σ Gil (F,Fi)BP : (4.5.11) 


m=0n=m p= 


Note that, while we may ultimately be able to express the scattered field using a small number 
of sources, this is a fundamentally different strategy from what has been applied heretofore in 
the simple independent dipole models. In the latter, each source responds only to the primary 
field striking it locally. However, here the p?””, i= 1,2,...,N;eq responding sources act together, 
not in response to local stimuli but to express the response of the entire object to the distributed 
excitation of the pmn mode. In addition, using the MAS-MAS/TSA numerical code the p?”"” can 
be generated and stored for any number of pmn spheroidal modes in the infinite series “45 11). 
Once this is done, then the truncation criterion can be determined from the input primary magnetic 
field easily. 
Overall, the entire GSEA approach can be described briefly as follows: 


mn 


1. For a given UXO amplitudes of responding magnetic dipoles p?”"", i = 1,2,...Nyeq rings 
are determined and sorted in the universal library for any number of basic spheroidal mode 
excitations pmn = 1,2,.... 


2. Once step | is done, then for a given sensor the primary field is decomposed into spheroidal 
modes, the spheroidal modal decomposition bpm coefficients are calculated, and the neces- 
sary number of spheroidal modes is determined. 


3. Use reduced set of sources (RSS) ρὲ to calculate EMI response for each pmn-th basic 


excitation. Scale each pmn-EMI response on the Dpmn coefficients and calculate the target’s 
complete response by just superposing responses of each basic excitation. 


The GSEA, which is based on the MAS-TSA and introduced here, can produce a target’s ultra- 
wideband frequency response. Thus, the proposed GSEA can be used directly to obtain an object’s 
TD EMI responses via convolution theorem without recalculating amplitudes of the responding 
sources. This makes the GSEA a unified model to treat both FD and TD data, and is attractive from 
a practical point of view, as many state-of-the-art EMI sensors (EM-63, EM-61, Zonge NanoTEM) 
are operating in TD. 

Let us briefly describe the important formulas required to compute the TD EMI response for 
a general current waveform /(t) flowing in a transmitter loop. By using the convolution theorem 
[29] the induced voltage in the receiver coil can be expressed as 


-- [Αἰ (t—t)I' (t) dt—A(O)I' (t) —A’(t)1(0), (4.5.12) 


where A(t) represents an object’s impulse response and the prime means the derivative with respect 
to time f. Equation (4.5.12) represents the TD response of an object to a general excitation current 
I(t)source. 
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4.5.5 Frequency-Dependent Magnetic Susceptibility 


There are three types of magnetic susceptibility that generate soils’ EM responses: (1) induced, 
(2) permanent/remanent, and (3) viscous remanent (VRM). The soil magnetic properties are deter- 
mined by the presence of iron and iron-oxide particles. Permanent remanent magnetization is the 
magnetization that exists in the absence of any applied field. Induced magnetization is the magneti- 
zation that arises in the presence of an external magnetic field, and viscous remanent magnetization 
(VRM) is a phenomenon that occurs when magnetization of an object placed in an external mag- 
netic field changes in a time relative to the applied field. This means that the object’s susceptibility 
is a complex frequency-dependent [45] 


X(@) = x'(@) + jx"(@) (4.5.13) 


where wgs the angular frequency, j is the unit complex number, and y'(@) and 7”(@) are the real 
and imaginary parts for the frequency-dependent magnetic susceptibility. There are various refer- 
ences [31 and references therein] that discuss many aspects of soil magnetic properties. The most 
common frequency-dependent complex magnetic susceptibility, assuming that magnetic relaxation 
time constants are uniformly distributed between times τὶ and 72, is modeled as follows [49]: 


1 | 


- 0 [1-- ral 
x0) = 10 (\— era Moet 


(4.5.14) 


where Yo is the DC value of the susceptibility. 
For magnetically susceptible ground in the present of a metallic object, the field that is mea- 
sured by the sensors contains two parts and can be written as 


zyines 


A” (@) =H" (w) +H! (o) (4.5.15) 


where Η" (ΟἹ) and Hq”! (@) are respectively magnetic fields produced by the magnetically sus- 


ceptible soil and the object. The magnetic field H ΟΡ] contains all interactions between the object 
and the susceptible host medium. For determining contribution of each magnetic field H*’ (@) 


and ΗΠ (@) in the measured field, a full EMI problem must be solved. In this section the soil is 
considered to be a uniform half-space, and interaction between the soil and the object is taken into 
account by using the image source method [49]. 


4.5.d Numerical Results 


This section presents some numerical and experimental data that demonstrate the GSEA as a uni- 
fied model for any FD or TD sensor configuration, and its applicability to complex data sets: novel 
waveforms, multi-axis, vector, or tensor magnetic, or electromagnetic induction data, or any com- 
bination of magnetic and EMI data; and ability to take into account the influence of conductive and 
magnetically susceptible geological soils on metal detectors. 


4.5.d.(1) The universal RSS ΑΒ it was discussed above, the reduced source set (RSS) p?”"” in Eq. 
(4.5.11) depends only on the target’s geometry and electromagnetic properties. To validate such a 
unique characteristic of the RSS, here comparisons between RSS modeled and experimental data 
are given for an actual UXO (81 mm) in both FD and TD. The data were collected by two EMI 
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systems that are widely used in the UXO discrimination community: (1) a wideband frequency- 
domain sensor (GEM3) developed by Geophex Ltd. [95] and (2) a time-domain instrument (EM63) 
developed by Geonics [96]. The data were collected on the U.S. Army Engineer Research and 
Development Center test stand site. 
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Figure 4.5.2: Frequency-domain EMI response for 81 mm UXO; (a) GEM-3 excitation; (b) nose up vertical, 
(c) nose up 45 inclination, (d) horizontal. 


Figure 4.5.2 shows the comparisons between GSEA and actual data in the frequency domain. 
In this case data were collected for the UXO oriented in three different directions relative to the 
GEM.-3 sensor’s head: (1) vertical tail up, (2) 45° degree nose up, and (3) transverse. The GEM-3 
frequency range is from 30 Hz to 50 kHz. The comparisons between measured and actual exper- 
imental data are in very good agreement for all orientations. Note that the RSS produces results 
in an ultra-wideband frequency range, from magneto-static (0 Hz) to EMI frequency limit. This 
allows users to accurately compute the scattered field at any required frequencies by simple inter- 
polation and to obtain EMI responses readily in TD via inverse Fourier transforms. To illustrate 
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this capability, here the first impulse responses in TD are evaluated by applying the digital filter 
technique [92] to the inverse sine transform [93] as follows: 


2 οο 
A(t) = =f ImB (a) sinot do (4.5.16) 
0 


where JmB (ΟἹ) represents an imaginary part of the magnetic flux that is calculated via RSS. Then 
the induced voltage is calculated using the time convolution technique (4.5.12). 
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Figure 4.5.3: Time-domain EMI response for 81-mm UXO; (a) ΕΜ-3 excitation; (Ὁ) nose up vertical, (c) 
nose up 45 inclination, (d) horizontal. 


In the EM63 instrument, the current waveform consists of an exponential current increase fol- 
lowed by a linear ramp off. The current has the three pulses per measurement. For comparisons 
between RSS and TD data, here the same 81-mm UXO is chosen. The object was excited in three 
orientations1: (1) vertical tail up, (2) 45° degree nose up, and (3) transverse orientations. For all 
three excitations, the vertical distance between the sensor’s transmitter loop and the center of the 
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cylinder is h = 60.00 cm. The TD induced voltage is calculated by an inverse Fourier transform of 
the frequency-domain magnetic field flux using equations (4.5.12) and (4.5.16). The comparisons 
are depicted in Fig. 4.5.3. They show very good convergence between the measured and modeled 
TD data. Therefore, the universal RSS should be able to produce EMI responses for a given target 
for any sensor in FD and TD domains. 
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Figure 4.5.4: EMI response for a magnetically susceptible soil versus the sensor height and frequency. (a ) 
Inphase part; (b) quadrature part. 


4.5.d.(2) EMI response for soil with frequency-dependent susceptibility In this section, first FD 
EMI response is studied for a magnetically susceptible half-space. The half-space is illuminated 
by a FD sensor. In these simulations for the sensor model the following parameters are used: 
current J, =1 A, and a 100-cm x 100-cm transmitter loop. The soil’s frequency-dependent sus- 
ceptibility is assumed to be the same as in equation (4.5.14) with a realistic χο = 0.005 D.C value 
of susceptibility [33], and τι = 10~° [sec] and τὸ = 10 9 [sec] time constants. Figure 4.5.4 shows 
soil’s responses in-phase (right) and quadrature (left) parts as a function of frequency and sensor 
height. The results illustrate that the soil’s EMI responses strongly depend on both the frequency 
and sensor height. As the sensor approaches the soil, soil’s response increases and it stays almost 
constant for antenna heights between 1 cm and 10 cm. Note that all parameters in Figure 4 are 
in Logarithmical scale. The soil response’s quadrature part approaches to maximum between 10 
kHz and 100 kHz, whereas at low frequencies the in-phase part is dominant. At highest frequencies 
(more than 100kHz) both parts of soil’s response decrease. Thus the soil with frequency-dependent 
magnetic susceptibility produces significant EMI responses almost entire UWB frequency range. 
Finally, to illustrate the soil effect on a buried object’s EMI responses, several numerical exper- 
iments were done in both FD and TD. The 81-mm UXO was buried under a magnetically suscep- 
tible half-space. The half-space electromagnetic parameters are exactly the same (%o = 0.005, τὶ 
= 10~° [sec] and τὸ = 10~? [sec]) as in the previous paragraph. The entire structure is illuminated 
by (1) the GEM-3 sensor with current J, =1 A, 10- and 20-cm radii coils, and (2) the EM-63 with 
I, =1 A, 100-cm x 100-cm transmitter loop. The sensors are placed 10 cm above the half-space. 
The UXO is oriented 45° nose up and its center is at 32-cm depth for the GEM-3 excitation and 
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50 cm for the EM-63 sensor. The results are depicted in Fig. 4.5.5. These results suggest that 
a magnetically susceptible half-space could significantly affect both frequency and time domain 
EMI responses. 
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Figure 4.5.5: EMI response for 81 mm UXO; (a) Frequency domain, (b) Time domain. 


These calculations are based on the functional form for ground susceptibility given by 4.5.14, 
which the Steve Billings, et al. at UBC/Sky state has “applied to every soil they’ve ever encoun- 
tered” in the TD. Also, the value of 5 x 10-3 is on the order of magnitude of the single frequency 
susceptibility values obtained by Janet Simms in numerous measurements at the UXO standardized 
test sites [33-35]. However the question arises whether this is the appropriate value to insert for 
%o. In particular, we have never seen FD soil effects of the magnitude and type shown in Fig. 4.5.5, 
left. This may mean that functional form 4.5.14 pertains successfully to TD, but that another larger 
effect appears in the frequency domain corresponding to the instantaneous soil response. The latter 
would not appear (or at least would only appear indirectly) in the TD measurements. The method 
used here to perform “absolute” calibrations of instrument data for conversion into actual magnetic 
field units allows us to infer actual broadband ¥ values from such data. Applying this to both the 
new GEM-3%* as well as emerging TD instruments should allow us to complete the characteriza- 
tion of soils of interest, modifying the functional form 4.5.14 as necessary to include all relevant 
FD and TD effects. 


4.5.e CONCLUSION 


The generalized standardized excitation approach, which is a fast, universal, and rigorous forward 
modeling system, has been developed and demonstrated. The GSEA is applicable to any FD or 
TD sensor configuration, and to any data set: novel waveforms, multi-axis, vector and tensor, or 
magnetic or electromagnetic induction data, or any combination of magnetic and EMI data. 

The proposed system has been tested against actual data in both the frequency (GEM-3) and 
time (EM-63) domains. Excellent agreements between the GSEA and experimental data have been 
demonstrated here. The GSEA technique could be used for building a new type of EMI sensor as 
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well as for optimal survey designing. 

By using an appropriate dyadic Green’s function the GSEA takes into account a host medium. 
EMI response from a magnetically susceptible half-space is analyzed versus sensor height and 
frequency. Numerical tests show that a frequency-dependent, magnetically susceptible half-space 
significantly affects the 81-mm UXO’s EMI responses in both FD and TD. This study suggests 
that a magnetically susceptible host medium must be taken into account in UXO discrimination 
problems. 


4.6 NSMS Modeling of GEM-39* and HAP Method 


The normalized surface magnetic source (NSMS) model [97-100] is a rigorous physically com- 
plete model of any metallic target’s interaction with an arbitrary primary field in the electromag- 
netic induction (EMI) regime. Based on a number of well-defined measurements for a given target, 
the system solves for its fundamental, characteristic response parameters. After target character- 
ization, the response of that target to any combination primary field modes can be quickly and 
accurately calculated. 

The NSMS model can be considered as a generalized surface dipole model [68, 82—87, 95, 96, 
101, 102], which in a special limiting case coincides with the infinitesimal dipole approach. In 
the model, the response of an object to the exciting field of a sensor is taken to originate from a 
layer of equivalent magnetic dipoles distributed over a spheroidal surface surrounding the target. 
Such a distribution can be generated by spreading positive charge density Q over the exterior of 
a thin smooth surface and an identical negative charge on its inner surface, resulting in a double 
layer of charge separated by an infinitesimal distance. This configuration introduces the proper 
discontinuities in the tangential components of the magnetic flux density vector B but does not 
affect its normal component. (A single layer of charge, on the other hand, would in no way affect 
the transition of the tangential components of the flux density but could model a discontinuity 
in the normal component of B.) Since there are no free magnetic charges in nature, the normal 
component of the magnetic flux density must always be continuous across the boundary between 
two media; the secondary magnetic field outside an object can thus be accounted for by a double 
magnetic charge layer of density Q distributed on an auxiliary surrounding surface; this density is 
proportional to the component of the primary (exciting) magnetic field normal to the surface. 

The amplitudes of the auxiliary sources are determined directly, and with great speed and accu- 
racy, by minimizing the difference between measured and modeled data for a known object-sensor 
configuration. The charge density Q is a property of the target, and its integral Q over the spheroidal 
surface constitutes a global magnetic “capacitance” of sorts that for a given object is invariant in 
the sense that different computational constructs — e.g., surrounding surface, measurement grid, 
object location or orientation — and different primary fields produce the same value. Once the 
source strengths are known, the NSMS model can be used in discrimination processing that either 
compares the measured fields to those stored in a library of known objects — 1.6., a “pattern match- 
ing” technique — or uses the integrated NSMS as a discriminant -- 1.6., “genuine” inversion. In 
“pattern matching” discrimination processing, an optimization or search algorithm is used to de- 
termine which known UXO has the catalogued ultrawideband Q distribution that best reproduces 
the signal received by the survey sensor. Used in this way, the NSMS system is a faster forward 
model than other physically complete forward models such as the standardized excitation approach 
(SEA) or the generalized SEA [88, 103-106]. 

We use the total NSMS, Q, of a target as a discriminant in Sec. 5. This Q may also be 
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parametrized by a Pasion Oldenburg type model and/or combined with other parameters, for ex- 
ample from the dipole model, in order to further characterize a target and classify it based on its 
EMI response [107]. When extracting Q from training data in order to create a given target’s “fin- 
gerprint,” the location of the target is generally know a priori. However, when the target’s position 
and orientation are not known, then inverse algorithm has to combine inverting not only for which 
target is present, but where that target is. This turns out to be a combination of a linear and non- 
linear process and is often time consuming. To this end, we have developed the HAP method (H 
(magnetic field), A (magnetic vector potential), and y (magnetic scalar potential)) which under 
certain assumption can find the location and orientation of a dipolar source without an expensive 
inversion [107-110]. 

While it’s easy to show that simple, infinitesimal dipoles are frequently deficient for modeling 
a target’s general response under diverse circumstances, the HAP approach assumes that under any 
particular circumstances a dipole representation is likely to furnish a good first order representation 
of the response at hand. In order to find an estimate of the target’s location and orientation in a 
rapid, robust manner, we temporarily assume that the sources of the anomalies being surveyed are 
caused by simple dipole sources. Using the magnetic dipole approximation in the EMI regime, the 
scattered magnetic field H, magnetic vector potential A, and magnetic scalar potential y can be 
written as: 


ot ee P), ae eee rete 


R2 47R3’ 47R3 


where R = 7—7’ with 7 and 7’ being the position of the observation and of the source, respectively. 
From (4.6.1), it can be shown also that y = (H-R)/2 and A = μι(Η x R). Now, taking R x A, 

we can solve for R as ΝΞ 

2yH - (A x H)/Ho 


Ap (4.6.2) 


Similarly, taking the cross product of A in (4.6.1) from the left side, the magnetic dipole moment 
P can be expressed as ᾿ 7 ΝΞ 
Ρ--4πκί(Εψ- (A x H)/wUo). (4.6.3) 


Thus, under magnetoquasistatic (MQS) assumptions of a dipole source, (4.6.2) and (4.6.3) repre- 
sent an analytical method for finding the location, R and orientation, P of the dipole. 

Currently, {H,A, yw} are found from EMI data through the Normalized Surface Magnetic Source 
Method (NSMS). An equivalent, flat, two dimensional NSMS surface located in between the mea- 
surement locations and the target can accurately reproduce the response from the target(s) being 
interrogated [100]. This NSMS surface consists of a layer of magnetic charges or dipoles from 
which {H,A, w} and hence {R,P} can be calculated. Note that no information regarding the tar- 
get(s) is required to calculate this NSMS surface layer. 

The NSMS model must be adapted to each unique instrument one time. This has been accom- 
plished for the GEM-39* and the results of matching GEM-39%* data to the predicted response 
from the NSMS model are shown in Figs. 4.6.1—4.6.3. 
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Figure 4.6.1: Comparison of H, from the GEM-3Y* and NSMS model predicted response. Data at 90Hz 


over a Steel sphere, inphase part. 


Benjamin Barrowes, Kevin O’ Neill -90- - ENGINEERING 
Ξ-- ΑΚ...) AT DARTMOUTH 


MM-1537, GEM-3Y*- Final Report 4 MATERIALS AND METHODS: GEM-39* SOFTWARE 
4.6 NSMS Modeling of GEM-3Y* and HAP Method 


'S0 (oa Predcted 
—GEM-3D plus data 


100 


Hy [PPM] 
Oo 


0 20 40 60 80 
Obser. Points # 


Figure 4.6.2: Comparison of H, from the GEM-39* and NSMS model predicted response. Data at 90Hz 


over a steel sphere, inphase part. 
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Figure 4.6.3: Comparison of H, from the GEM-3%* and NSMS model predicted response. Data at 90Hz 
over a steel sphere, inphase part. 3 levels with 81 points on each level. 
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4.7 Electromagnetic induction from highly permeable and conductive ellipsoids under ar- 
bitrary excitation — Application to the detection of unexploded ordnances 


Symmetrical objects such as spheres or spheroids may be useful as models of some subsets of 
possible targets; and for testing instruments they may be particularly useful when one can examine 
signals for expected symmetries. However, particularly when all vector components are being 
recorded and examined, a truly 3-D object solution offers the most general reference response. 
Heretofore there has been no broadband 3-D analytical solution available, except by relatively 
laborious numerical methods. The approximate solution produced here is rapidly evaluated under 
diverse conditions. It represents a milestone in benchmark cases that can be used hereafter in a wide 
variety of contexts. Specifically in our discrimination processing, with this new solution we can for 
the first time model the responses of unseen objects quickly without having to assume restrictive 
symmetries. We now have a 3-D solution to guide analysis, statistical learning machines, etc, for 
inferring overall geometrical properties such as volume and aspect ratios, from signals across the 
breadth of the EMI band. 

The secondary field produced by three-dimensional highly permeable and conductive objects 
is computed in the electromagnetic induction regime, with the purpose of modeling unexploded 
ordnances and surrounding clutter. The analytical formulation is based on the ellipsoidal coordi- 
nate system, able to model real three dimensional geometries as opposed to bodies of revolutions 
like within a spheroidal approach. At the frequencies of interest (tens of Hertz to hundreds of 
kilo-Hertz), the conduction currents in the soil are negligible and the fields are computed in the 
magneto-quasistatic regime based on the Laplace equation. Inside the objects, where the wave 
equation governs the field distribution, the currents are assumed to have a small penetration depth, 
allowing for the analytical simplification of the field components, which become decoupled at the 
surface. This approximation, valid across the entire frequency spectrum because of the high per- 
meability and conductivity, avoids the necessity of using ellipsoidal wave functions and results in 
a considerable saving of computational time. Numerical results compare favorably with numerical 
and experimental data . Finally, the optimization approach used to match our numerical predictions 
with experimental data demonstrates the possibility of remotely inferring the material properties 
of objects. 


4.7.a Introduction 


One theoretical approach models the UXO as a (possibly) magnetic and (possibly) conductive 
spheroid [2, 111-113]. The model is therefore mathematically more involved but offers impor- 
tant flexibility in terms of geometries (from elongated to flat objects within the same formalism), 
material parameters, and distance to sensor. 


In this section, we extend the second approach to ellipsoidal geometries with the underlying 
motivation of providing an improved physical model of the clutter surrounding UXO, directly use- 
ful for example in statistical and optimization methods such as in [114-118]. Clutter is indeed 
most often composed of random three dimensional (3D) pieces of steel which are jamming the 
measured signals and are in part responsible for the high false positive alerts. Their proper mod- 
eling is therefore essential and the ellipsoidal coordinate system provides a general framework to 
do so, allowing for the modeling of real 3D geometries while retaining the speed advantage of 
analytical calculations, fundamental for example in inversion algorithms [1 19][120][121—123]. In 
addition, we use in this paper the small penetration approximation (SPA), justified by the major 
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steel component of clutter objects. SPA was originally developed as a high frequency approxima- 
tion for spheroidal geometries [113] under the assumption that the currents induced on the surface 
of the object penetrate only slightly into its volume. Remarkably, however, SPA has been shown 
to yield accurate results over the entire EMI frequency range for highly permeable and highly con- 
ductive objects, such as those made of steel. This success has prompted various implementations 
of the SPA, both analytical [113] and numerical [124, 125]. In this paper, we generalize the SPA 
to ellipsoidal geometries and confirm the good agreement obtained over the entire frequency range 
by comparing our predictions with both well established numerical results and with measurement 
data. 


The configuration under study is depicted in Fig. 4.7.1: a homogeneous ellipsoid of semi-axes 
(a > b >c) along the directions (£,9,2) (see the Appendix for mathematical details), of perme- 
ability μι = 100uUp9 and conductivity σὶ = 10’ [S/m] is embedded in a non-magnetic and non- 
conductive soil which, at the frequencies of interest, is essentially transparent and modeled as 
free-space. A sensor located above the object emits a primary field Η΄ and measures the secondary 
field H’ produced by the object. We subsequently use two models for the primary field: a sim- 
ple uniform field, valid for small or deeply buried objects, and the field produced by the GEM-3 
instrument [1], more realistic for large or close objects. 

In the EMI regime, the magnetic field outside the ellipsoid is irrotational and can therefore be 
written in terms of a potential function which obeys the Laplace equation. Denoting the primary 
and secondary potentials as U° and {75, respectively, the total magnetic field in the region exterior 
to the ellipsoid is written as 


H. =H’ +H =-VU°-VUs. (4.7.1) 
In the interior region to the object, the magnetic field H; obeys the wave equation 
VxVx A, -KH, =0 (4.7.2) 


where ki Ξξειῶσι μι. This limiting form of k; that applies in this parameter space means that (4.7.2) 
leads in fact to a diffusion equation representing magnetic diffusion, though it retains the form of 
a wave equation. The two quantities H; and H> are related by the usual boundary conditions. 

For details on the ellipsoidal coordinate system, the separation of the wave equation in ellip- 
soidal coordinates, and the derivation of the secondary field H” in the ellipsoidal coordinate system 
as a function of the type of incidence (uniform and GEM-3) within the SPA approximation, see 
[126]. 

The ellipsoidal coordinate system is defined based on the definition of the master ellipsoid of 
semi-axes dimensions (a,b,c) along the (£,,2) directions respectively (a > b > c): 

2 ak go 
xy ζς 
szt+a+5=1. de 
Upon defining the variable k? = a* — c? and h? = a* — b* [127], three families of orthogonal sur- 
faces can be obtained [128], used to defined the ellipsoidal coordinate system (p, U,V): 
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Figure 4.7.1: 


Configuration of the problem: an ellipsoid of material parameters (μι = 100U0, o1 = 


107 [S/m]) is surrounded by free-space (2 = Uo, O2 © 0) and is excited by a GEM-3 instrument [1] modeled 
as two current loops (figure not to scale). The semi-axes of the ellipsoid are denoted by a, b, and c (with 
a>b>c) along the ἃ, ¥ and 2 directions, respectively (see the Appendix for mathematical details on the 


ellipsoidal coordinate system). 


and where p? > k? > pu? >h? >v?>0. 
The magnetic field is obtained from the gradient 


1 dU 1 οὔ 1 οὔ 
VU=) (i v 4.7.5 
Php Op ὦ Ὁ ομ “hy ὃν ΡΝ, 
and is written as 
H’ =-VU° =-p > ΣῈ (p’, u,v) — τ ΣΡΡῈ (ρ, μ'. νὴ -- ὩΣ E(p,y,v’) 
hp n,p 

(4.7.6) 

for the primary field. 
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Introducing the field expansions and separating the terms functions of ΡΨ ) and ΒΡ) yields 


(p) J H2 py ͵ hp O (hy , 
LBr {i [1 (a)E(a) +1(a)E'(a)] E(u, νὴ — ii on (Ft) I(a)E(a,p',V) 


h hy 9 {ἢ 
Ε΄ 305 ra (1) Τ(ΛΕ(α, μ,ν) — τ: ful αΛΕία, μεν ᾿ 
Ξ (p) μ / hp 0 ( ) ! hp " 
=) δὴ — ik; —E(p’,U,Vv + —— | — E(a,w,v)+—fyE(a,u',Vv 
Σ [ ἱκι---Ε(ρ΄,μ " ἬΝ hy lt (α, μ΄, ν) 27} (α, μ΄ ,.ν) 
hp 0 ͵ hp " 
τ ὅν (Fete) Elen) + νίαν" 417) 


and it is understood that this expression is evaluated at the surface of the ellipsoid p = a. It should 
be noted that most of the terms in Eq. (4.7.7) are straightforward to compute: the derivatives and 
second derivatives of the Lamé functions with respect to μ᾽ and v do not present any extreme 
difficulty. 


4.7.b Secondary fields due to uniform incidence 


The first situation we shall consider is the one of a uniform incidence. The latter is not merely a 
theoretical idealization since objects that are far from the GEM-3 sensor, or those that are small 
(in the order of a few centimeters), do not see a very inhomogeneous primary excitation. The 
assumption of uniform incidence in these cases is therefore excellent. 

For the sake of later comparison with known solutions, we first suppose the primary field to be 
in the (xz) plane, although the general case can be treated just as easily. A uniform field derives 
from a simple linear potential in x and z, which can also be expanded in the ellipsoidal coordinate 
system as 

U° = Hoxsin€ + Aozsind = Σ ΔῈ P) (7) (9 μιν) (4.7.8) 
np 
where Θ = 0 for axial excitation and θ = 7/2 for transverse excitation, and where Ho is the am- 
plitude of the uniform field. Realizing that x and z can be directly expressed in terms of Lamé 
functions of the first kind with n = 1 we write 


Ao sin@ __(1) sn OOO cos = Eo 
= Ε; 4.7.9 
so that the bl? ) coefficients of Eq. (4.7.7) are directly given by 
a) _ Hosin@ (2) (3) Ho cos 0 (p) 
b= b,’ =0, by’ =——=—, dD’, = 0. 4.7.10 
1 hk’ 1 ’ 1 κν-. "3 n#l ( ) 


The expansion of the secondary field also deserves attention before proceeding with the numer- 


ical solution of the ΒΨ ) Indeed, there exist an analytical solution to the external potential due to a 
homogeneous ellipsoid [129]: 


2 2 
y Ζ ds 
V (x,y,z) =f -- "ΕΝ πον FIFE πε΄ (4.7.11) 
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where ko is a constant. We see that the first term of Eq. (4.7.11) is directly proportional to Ko (i.e. 
n =), while the successive three terms are directly proportional to Κι, 1, and M, (i.e. n= 1), 
respectively. The potential of Eq. (4.7.11) can therefore be expressed with exactly, and only, four 
ellipsoidal term, directly related to the demagnetization factors discussed in [130, 131]. In addition, 


since the n = 0 term is discarded in our solution (the ellipsoid being uncharged), the expansion of 


the secondary potential is reduced to n = 1 only (p = 1,2,3), and Be, = 0. Moreover, in the 


specific case of a uniform primary field in the (xz) plane, B°) = 0 and we are left with only two 


unknowns to determine. 


We first apply our method to the computation of the secondary field produced by a series of 
spheroidal geometries, under axial, transverse, and oblique incidences. The main purpose here is to 
being able to compare the results with those obtained in our previous work [2]. The correspondence 
between the elongation of the spheroids and the size of the ellipsoids is given in Table 2. The axial 


Spheroidal Ellipsoidal dimensions Measured field 
elongation a>b>c (along &, ¥, 2) Fig. 4.7.2(top) | Fig. 4.7.2 (bottom) 
e=6 a=1.5cm, b=c=0.25 cm A, H, 
e=3 a= 1.5 cm, b=c=0.5cm A, H, 
e= a=b=c=1.5cm Η, or Hz 
e=1/3 | @a=b=15 cm, c=0.5 cm Hy Ay 
@=1/6 | a=b=15 cm, e= 0.25 cm Hy H, 


Table 2: Correspondence between the elongations of spheroids and the size of the ellipsoids as shown in 
Fig. 4.7.1 and defined by Eq. (4.7.3). 


and transverse cases are shown in Fig. 4.7.2 for spheroids located 10 cm below the sensor and 
with elongations e € {6,3,1,1/3,1/6} (e > 1 corresponding to prolate objects while e < 1 to 
oblate ones). The overall responses of the objects as function of frequency follow the now well- 
understood trend: at low frequencies the secondary magnetic field is aligned with the exciting field 
while at high frequencies, the induced currents are confined to the surface of the ellipsoid and are 
out of phase by 7 due to Lenz’ law. The imaginary part of the field therefore vanishes at both ends 
of the frequency spectrum and peaks at an intermediate frequency corresponding to the resonance 
of the object. 

In the axial case, the magnetic field sees more object surface with prolate geometries, produc- 
ing a stronger response than oblate geometries, while this tendency is reversed under transverse 
excitation for the same reason. The position of the peak in the imaginary part is directly related 
to the resonance of the object, also closely related to the size of the object along the direction of 
the primary field. In the transverse case, elongations of e € {1/6,1/3,1} present the same size 
along the primary field and therefore, a similar location of the resonant peak. Prolate objects, on 
the other hand, have a varying size along this direction and thus, different resonant peaks. Similar 
arguments can be applied to explain the results obtained under the axial excitation. 


The first conclusion when examining Fig. 4.7.2 is that the agreement over the entire frequency 
range is reasonably good in most of the cases, thus proving that the SPA is an adequate approxima- 
tion for highly magnetic and conductive objects, as already shown in [113]. The prolate cases are 
seen to usually agree better with the reference markers [2], and even more so under the transverse 
excitation. The worst agreement is obtained for oblate geometries under transverse excitation, with 
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increased disagreement as the objects flatten, in agreement with the conclusions reported in [125]. 
Although not shown for the sake of clarity, simulations have been extended up to extreme elon- 
gations of 10 and 1/10 and confirm these findings. A possible explanation is that under these 
circumstances, a very large portion of the object surface has a very high radius of curvature so that 
a dominant radial decay approximation of the currents may become less accurate. 


Fig. 4.7.3 presents the results for a configuration similar to the one corresponding to Fig. 4.7.2 
but where the object is rotated at 45 degrees. Similar conclusions to the previous case can be 
drawn: prolate spheroids are well modeled whereas oblate ones present the highest mismatch, the 
latter being more pronounced for more extreme geometries. The disagreement for e € {1/3,1/6} 
is here seen to be more evident at low frequencies, while at high frequencies (where the SPA is 
more and more valid) the agreement is usually reasonably good. 

Finally, Fig. 4.7.4 shows the secondary field produced by a real ellipsoid of dimensions 1 cmx2cmx3 cm 
along the x, ¥, and Z directions, respectively, under an axial excitation. The depth and material pa- 
rameters are otherwise identical to the previous cases. The reference solution is obtained from 
the method of auxiliary sources (MAS) presented in [132],[3], and has been truncated at 1 MHz 
because of the fine mesh and long computation time required at higher frequencies. At low fre- 
quencies, however, the MAS is very accurate and shows that our method suffers from some errors 
which we attribute to the choice of the location of the matching points. Nonetheless, the difference 
between the secondary field produced by an object of dimensions 1 cmx2 cmx3 cm (the present 
ellipsoid) and 1 cmx 1 cm x3 cm (a spheroid of elongation e = 3) is quite significant and allows 
us to conclude that away from the lowest frequencies, the agreement between the two methods is 
reasonably good. 


4.7.5 Secondary fields due to GEM-3 incidence 


Although the assumption of uniform primary field is sufficient in many cases, it is still beneficial to 
model the primary field generated in real field measurements, typically as obtained from the GEM- 
3 instrument. The GEM-3 is a broadband electromagnetic sensor manufactured by Geophex [1], 
whose transmitting sensor head consists of two current loops of radii of 20 cm and 11.074 cm. 
These two current loops are arranged as bucking coils with the smaller loop placed concentrically 
in the interior of the larger loop (schematically illustrated in Fig. 4.7.1). The two currents are 
optimized so that the primary magnetic field at the center of the instrument is minimal and the 
measurements closely correspond to the secondary field ΗΝ. 

The analytical expression of the axial and transverse magnetic field components produced by a 
circular current loop in a cylindrical coordinate system, based on the complete elliptical integrals 
of the first and second kind, can be found in [74, 75, 133] and reads: 


Holo (z—2z0) γ2-ὐ2- (ς-- 29)Σ 
Βρ(γ, φ,1) = "στ τ ΡΣ K(ke) + Gay ae pete) (4.7.12a) 
= Lolo r —re+(z—Z0)" 
Blt .2) = a κι ΥΠΠΠ᾿ IK (ke) - ark GaP tte): (4.7.12b) 


where /o is the radius of the current loop located in the z = zo plane, J, is the current flowing in 
the @ direction, and 


4ar 
ko= lar Cas (4.7.13) 
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Figure 4.7.2: Real and imaginary parts of the secondary field H” produced by spheroids of elongations 
shown in the legend, buried 10 cm below the sensor, and under axial (upper panels) and transverse (lower 
panels) excitations. The material parameters are μι = 1009 and σι = 107 [S/m]. The sphere case (denoted 
by ’e=1 (n)’) has been normalized by 1/3 for display convenience. The markers are obtained from the 


method presented in [2]. 
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Figure 4.7.3: Similar to Fig. 4.7.2 for an object buried at 30 cm and rotated by 45 degrees. Upper panel: Hy. 
Lower panel: H,. Note that for this specific case only, the ranges have been set to v € [0,0.7]h for prolate 


spheroids and yt € [0,0.7](k —h) +A for oblate spheroids. 
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Figure 4.7.4: Real and imaginary parts of the secondary field produced by an ellipsoid of dimensions 1 cm 
in ἃ, 2 cm in ¥, and 3 cm in Z. The material parameters are μι = 100Up and σὶ = 107 [S/m], and the ellipsoid 
is buried 10 cm below the sensor. The reference markers have been obtained using the method of auxiliary 
sources [3]. 


K and E are the complete elliptical integral functions of the first and second kind, respectively. 
Note that this model agrees with the alternative methods of using a numerical integration or the 
application of the Biot-Savart’s Law, and has the added benefit of being computationally less ex- 
pensive. 


The first step in modeling a GEM -3 incidence is to determine the expansion coefficients ΡΨ ) of 
the primary field in Eq. (4.7.6), simply done using a point matching technique. Hence an ellipsoid 
is first defined in space, yielding a series of (x,y,z) components. At those locations, the magnetic 
field is computed from Eq. (4.7.12), and the field components are transformed into (Hp, Hy, Hy). 
Eq. (4.7.6) is then used across many points to compute the matrix of an over-determined system, 


from which the pl? ) are solved for in the least square sense. In doing so, we find that the coefficients 
are more stable when the matching is done on Hy for prolate-type objects and on H,, for oblate-type 
objects. The matching using Hp, reasonably good in both cases, can also be used for the purpose 
of comparison. 


Due to the non-uniform nature of the primary field, terms of higher degrees than n = 1 need to 


be included in the ellipsoidal expansion and therefore, all the BY, coefficients need to be solved 


for simultaneously. The lack of direct orthogonality between the various terms of Eq. (4.7.7) results 


in a somehow unstable numerical solution: the influence of some modes (n, p) can leak into other 


modes and degrade the accuracy of the corresponding BY ). This is typically the case when some 


BW?) are very small, and need to be solved as zero numerically: although the results are numerically 


small as expected, the non-zero BY ) can be slightly modified, introducing errors in the secondary 


field values. In the previous examples, this phenomenon has been bypassed by using a physical 
argument to force some coefficients to zero and to solve for the non-zero ones only. In the non- 


(p) 


uniform incidence case, however, such reasoning cannot in general be applied and all the B,; ’ need 


to be solved for simultaneously. 


In order to minimize the numerical leakage problem, we resort to solving for the BY ) in two 


steps. First we solve for all the BY? ) in a usual manner. Second, the ΒΡ ) are used to compute 


the field at the sensor location. The contribution of each mode is compared and those modes 
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Lod, g(trequency) fz] Lod, q(trequency) fz] 
Figure 4.7.5: Real and imaginary parts of the secondary field produced by a spheroid of elongation e = 3 
(left panel) and ὁ = 1/3 (right panel) buried 10 cm below the GEM-3 instrument. The material parameters 
are μι = 1000 and σι = 10’ [S/m]. The primary field is modeled using Eqs. (4.7.12). The reference 
markers (marked ‘ref’ in the legend) are obtained from the method presented in [2]. 


contributing less than a threshold (1/1000 of the maximum contribution in our case) are discarded. 


The BY ) are then solved anew only for the contributing modes. Note that this selection works best 
if it is performed at the first frequency, and kept identical for subsequent frequencies. 


The results for a GEM-3 axial excitation of a prolate (e = 3) and an oblate (e = 1/3) spheroid 
are shown in Fig. 4.7.5. The necessity of including more modes in the ellipsoidal expansion is 
apparent, although the number of modes is still significantly lower than under a spheroidal coordi- 
nate approach. The agreement obtained here with up to n = 3 for the prolate case and up ton = 2 
for the oblate case is seen to be good. The transverse excitation exhibits a similar good agreement 
and is not shown for the sake of brevity. 


Finally, measurement data are taken as a last validating element of our analytical approach. 
Measurements were performed using the GEM-3 instrument along the three axes of a steel ellipsoid 
of size 15 cmx7.5 cmx3.75 cm in the *, ¥, and Z directions, respectively (yielding an aspect ratio 
of 1:2:4), The magnetic permeability and electric conductivity of the steel are a priori unknown 
and can be optimized for when matching the experimental data to the theoretical ones. 

The comparisons between measurements and analytical results are illustrated in Fig. 4.7.6 for 
the data collected along the three axes. The caption provides information on the scaling factors y 
optimized when matching the data as well as on the values of the permeability and the conductivity 
for each case. It can be seen that overall, the numerical approach is very successful at predicting 
the EMI response of the ellipsoidal object across the entire frequency range. This good agreement 
represents an important validation of the method presented in this paper. 

A more specific examination of the results reveals both good and less good aspects, however. 
First, the scaling factor, which converts the units of GEM-3 measurements into [A/m], is relatively 
constant in the three independent optimizations results and yields y ~ 4.5 x 10°. Along with the 
overall good agreement in the magnetic field, this is an important result which confirms the ac- 
curacy of our measurement technique. Second, the first and second cases in Fig. 4.7.6 produce 
very similar results for the permeability and the conductivity, which is another important positive 
result. The numerical values confirm the expectation that steel is a highly permeable and conduc- 
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Figure 4.7.6: Simulated (solid lines) and experimentally measured (circles) real and imaginary parts of the 
magnetic field along the three axes of an ellipsoid: case A is along the longest dimension (15 cm), case B is 
along the second longest dimension (7.5 cm), case C is along the shortest dimension (3.75 cm). The location 
of the resonant frequency in the three cases follow this geometrical progression, as expected. The optimized 
parameters in the respective cases are: (μι = 140p9, σὶ = 4.5 x 10° [S/m], y = 4.7 x 10%) (μι = 1509, 
σι =5 x 10° [S/m], y= 4.5 x 103) (μι = 209, σι = 1.2 x 10° [S/m], y = 4.2 x 103). The dashed line in 
case C has been obtained with (μι = 100Up, σὶ = 1 x 107 [S/m], y=42~x 10°). The numerical data have 
been obtained with tree modes in all cases (N = 3). 
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tive material, which is a necessary condition for our SPA approximation to be trustable across the 
entire frequency range. More importantly, these results show that we are able, to some extent, to 
infer the material properties of an object in a totally remote and non-invasive way. Third, it is seen 
that the last situation in Fig. 4.7.6 stands out from the other two, with different estimates of the 
permeability and of the conductivity. We believe that this case is tantamount to cases of very flat- 
tened oblate spheroids studied previously (e = 1/6), where the currents flowing around the object 
see a very large radius of curvature at all points around the circumference, producing less accurate 
results within an SPA approximation. Moreover, the cost-function associated with this third situa- 
tion presents many local minima so that the optimized values μι = 20Up and o; = 1.2 x 10° [S/m] 
have to be understood as being one possible solution. In fact, the dashed lines correspond to the 
values of the magnetic field obtained with μι = 100uUp and o; = 1 x 107 [S/m] (the initial guess 
in the optimization procedure), and yield a reasonably good match with the measured data as well. 
Consequently, we can consider this last case as less reliable than the other two, but yet constitutes 
an important situation to consider in future work. 


4.7.d Conclusion 


The motivation of modeling clutter surrounding unexploded ordnances (UXOs) as three-dimensional 
steel objects in the electromagnetic induction regime has prompted the study of the ellipsoidal co- 
ordinate system presented in this paper. The hybrid approach, combining a static formulation based 
on the Laplace equation in the region exterior to the object and a full-wave approach based on the 
wave equation in the region interior to the object, has been formulated based on the Lamé functions 
and the small penetration approximation (SPA). The latter, originally derived as a high frequency 
approximation, is shown to be valid across the entire frequency spectrum of the GEM-3 instru- 
ment for highly permeable and highly conductive objects, in agreement with previous conclusions. 
Importantly, the SPA bypasses the necessity of using ellipsoidal wave functions inside the object 
and results in a considerable reduction of computational time, yielding an algorithm that produces 
broadband secondary magnetic field responses in seconds under uniform excitation, and in tens of 
seconds under GEM-3 incidences. 

Numerical results have been favorably compared with three independent results: those obtained 
with a well-established method for the study of spheroidal geometries (treated here as a special 
case of an ellipsoidal geometry), those obtained with the method of auxiliary sources for metallic 
objects (where the ellipsoidal geometry is a particular case of a general three-dimensional shape), 
and those obtained by GEM-3 measurements of a manufactured ellipsoid. The good agreement of 
all these results prompts us to conclude that the method presented in this paper is both fast and 
accurate, and can therefore be used in subsequent inversion algorithms for the classification of 
buried UXOs in clutter contaminated soils as well as for the estimation of material parameters in a 
non-invasive way. 
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5 Results and Discussion -GEM-39* 


The GEM-39%* can be used in several different measurement modes, namely, dynamic mode, 
static mode, and hybrid mode. These correspond to what type of information the user wants to 
know. Dynamic mode may be used mostly for the detection of metallic objects, while static mode 
could give higher fidelity lower noise data for discrimination. The hybrid mode, which allows for 
measurements taken while the sensor head is in motion, and relies on post processing to extract 
individual data points, may turn out to be the most commonly used mode. 

In static mode, the GEM-3Y* would be placed on a certain location and held immobile while 
one data shot is taken. Because the sensor is not moving, the data may be integrated for a very 
accurate low noise data point. Several of these in a rough grid would aid inversion efforts toward 
discrimination. The main drawback to this mode is that the data acquisition sequence would require 
a longer time than data acquired in a dynamic or hybrid mode. 

Acquiring data in the the dynamic mode would allow the user to hold and rotate the sensor 
head on the handle much more freely without having to pause between the data shots. At a 10 Hz 
measurement rate, substantial amounts of data could be acquired in a short period of time. One 
concern when acquiring data in dynamic mode is the effects of sensor head movement on data 
quality. Another concern would be data stability of each 10 Hz measurement without integrating 
to reduce noise. If the sensor head is not moved too quickly, movement issues can be kept to a 
minimum. Furthermore, when a post processing step is added (resulting in measuring in the hy- 
brid mode) data acquired while the sensor head was moving too rapidly could be flagged and or 
discarded. To address stability concerns of the 10 Hz measurements, we plotted the standard devi- 
ation of a few seconds of data acquired while the instrument was not moving (see Fig. 5.1.11). The 
GEM-39* has exceptional stability over short time periods even at the 10 Hz measurement rate 
and especially in the mid frequency band. Note that we normally only regard signals greater than 
100 ppm as significant. During post processing, measurements taken sufficiently close together 
(for example, when the instrument was held immobile during dynamic acquisition mode) could be 
averaged together leading to higher quality, higher SNR data. All of these data collection modes 
rely explicitly on a very accurate positioning system such as the beacon system described above 


5.1 Initial Data and Reasonableness Checks 


We have taken some initial data shots of a spherical object. The main purpose of this initial data set 
was for a reasonableness check. That is, do the separate components of the received secondary field 
have the right sign, reasonable amplitude, range, etc. Figure 5.1.1 shows the path of a sphere as it 
is moved along a line below the sensor head and is offset in the positive y direction. Figure 5.1.2 
shows the result if of extracting the quadrature data recorded at 450Hz and plot it as the sphere is 
moved along this path. Reasonableness checks on this data include 


e H, is always negative 
e H, starts negative then switches to positive 
ὁ H, is always positive 


e |H,| > |H;| due to the proximity of the sphere to the x — y plane 


results are symmetrical around the yz plane 
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Figure 5.1.1: Initial reasonableness check of the GEM-3Y*. Object is moved along a path below and offset 
to the GEM-39*. 


Figures 5.1.3 and 5.1.4 show shots from along a line crossing the center of the 3 GEM-39* 
Rx coils, on a 45° angle relative to the Hx and Hy coils. Spectra are shown for three positions: 
the one (10 cm) before crossing the origin (’pre’’), over the center, and after (’post’’) crossing the 
origin. Symmetry & appropriate pattern (magnitude) similar to those noted above are evident. 

Figures 5.1.5 and 5.1.6 show the components from Rx coils 1 & 2 (horizontal: Hx and Hy). 
Values are pretty small, not far above the noise/background variation floor (the sphere was quite 
small). The responses pass both quantitative and quantitative reasonableness checks based on the 
amplitudes and signs of the 3 components. 
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Figure 5.1.2: Initial reasonableness check of the GEM-3Y*. Quadrature response at 450Hz. 


Z component at center 
100 7 


50 - PF tea chia 

Inphase 
& ᾿ 

Quadrature ῃ - inate ca 

(ppm) 


Inphase 


-200 ἢ 
100 1000 104 108 


Frequency (ΗΖ) 


Figure 5.1.3: GEM-3Y* Z component of received field with sphere directly over the center of the instru- 
ment. 
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Figure 5.1.4: GEM-3Y* Z component of received field with sphere to the sides of the instrument. 
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Figure 5.1.5: GEM-39* X and Y component of received field at 3 positions. Inphase. 
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Figure 5.1.6: GEM-3Y* X and Y component of received field at 3 positions. Quadrature. 
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We also acquired data over 3in spheres of brass, aluminum, and steel in a 9x9 grid. Figures 5.1.7 
and 5.1.8 show data from different transects of the received field for the H, and H, components 
respectively for the steel sphere. These data traces should be on top of each other if the sphere was 
exactly aligned. Since the data are slightly off but in a consistent direction, we conclude that we 
are in fact seeing the results of a slight misalignment from center. 
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Figure 5.1.7: H, component of data from different transects of the received field over a 9x9 grid. 


Figures 5.1.9 and 5.1.10 show the H, and the H,, Hy components of the received magnetic 
field respectively over a steel sphere located southeast of the sensor. As expected the transverse 
components of the received magnetic field are negatives of each other. Also, all three components 
of the received magnetic field are available across the entire band from 30 Hz 100 kHz. 

Figure 5.1.11 shows the standard deviation of the measured H from the GEM-39* combining 
10 seconds of 10Hz data. The values are seen to be quite stable over the short time periods of 
collecting one data show. The GEM-39* does exhibit some slight drift during the warm up phase 
(15 minutes) especially, similar to other GEM’s. 
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Figure 5.1.9: H, component of the received magnetic field over a steel sphere located the southeast of the 


sensor. 
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Figure 5.1.10: H,, Hy components of the received magnetic field over a steel sphere located the southeast of 
the sensor. 


(ppm) 


std(H) 


0 | | | 
10° 10 
frequency (Hz) 


10 


Figure 5.1.11: Standard deviation of the raw 10Hz data of the 3 field components of the received magnetic 
field. 
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Figure 5.2.1: GEM-39* on the lab grid for blind test suite #1. 


5.2 Blind Test Suite #1 Training Data 


Blind tests were performed for the GEM-39" in that laboratory data was acquired at CRREL and 
handed to our colleague Fridon Shubitidze at Dartmouth College. There are two stages to blind 
testing with the NSMS or SEA models. The first stage involves acquiring characterization data of 
a target in order to set up a library of responses from possible targets in the test. The second stage 
is inverting for position, orientation, and which target is in each test. These were single target tests. 

The characterization phase for this round of blind tests was complete in April, 2008. Data 
was taken on a 3x3 (20cm spacing, blue dots in Fig. 5.2.2 for background data) or a 9x9 (Scm 
spacing) grid at multiple elevations and orientations of each library target. Figure 5.2.1 shows 
the GEM-3%~ on the lab measurement grid, and Fig. 5.2.2 shows the grid numbering system 
used for this data. Figure 5.2.3 shows a 60mm mortar in the 3 test configurations used to acquire 
characterization data. Generally, the three 9x9 measurement grids were taken at: 


e horizontal, nose north (or +y), abbrev HNN. 
e nose north 45degrees down from horizontal, abbrev NND 
e nose north 45degrees up from horizontal, abbrev NNU 


Depth measurements are to the nearest point with the center of the target centered at (0,0), and 
background grids of 3x3 were acquired in between, before, and after each 9x9 measurement grid. 
The suite of targets we acquired characterization data for this test are pictured in Fig. 5.2.4. 
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Figure 5.2.2: Grid layout used for blind test suite #1. 


Figure 5.2.3: 60mm mortar in measurement setup for blind test suite #1. 
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Figure 5.2.4: Possible targets in library for test suite #1. 


5.3 Blind Test Suite #1 Results 


A library was constructed based upon the lab characterization data described above. The the HAP 
method and the NSMS model (see Sec. 4.6) were applied to subsequent blind test data (also on 
a 9x9 grid). The total normalized magnetic charge, Q for each blind test target was compared 
to the library Q and ranked according to best match. Figures 5.3.1—5.3.8 show Q for the blind 
test data compared to that of the best match from the library. In summary, Fig. 5.3.9 shows the 
geometry and tabulates the results from this blind test. All 8 targets were identified successfully 
with positions and orientations very close to the actual. This blind test is encouraging in that it 
suggests the NSMS model, combined with the HAP method for fast target pinpointing, is very 
effective at correctly discriminating between targets with vector data such as that provided by the 
GEM-39". 
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15000 


= ©- Quadrature Blind test #1 
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Figure 5.3.1: Total NSMS of blind test 1. 
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Figure 5.3.2: Total NSMS of blind test 2. 
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Total NSMC [Arb] 


Total NSMC [Arb] 
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Figure 5.3.3: Total NSMS of blind test 3. 
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Figure 5.3.4: Total NSMS of blind test 4. 
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Figure 5.3.8: Total NSMS of blind test 8. 
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Measurement surface 


Figure 5.3.9: Geometry diagram and positions for blind tests derived from the HAP method compared to 
truth. 
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November 29, 2007 


Figure 5.4.1: Contents of test plot #2 blind test suite #2. 


5.4 Blind Test Suite #2 -CRREL Test Plots with Beacon Positioning 


On October 23, 2008, CRREL personnel acquired data over 15 target clusters in the CRREL test 
plots. Each target cluster consisted of either one or 2 emplaced targets. The test plot emplacements 
consisted of 5 target clusters in 3 separate 20’x20’x8’ concrete plots. The items in each cluster 
are shown in Figs. 5.4.1—5.4.3. In these figures, a red box around an alternate target indicates a 
target that used to be in that position, but was thereafter replaced with the target listed. Target 
designations such as LG-22 and so forth are CRREL designations and refer to different pieces of 
scrap with specifications available upon request. Test plot #1 was not used for this suite of tests. 

The GEM-3%~ including the functioning beacon positioning system was used to acquire 
data over all 15 target clusters. Appendix C.3 describe the data processing of this data while 
in Sec. 5.4.a, the results of a multi-object NSMS inversion [98, 99, 134] is presented. References 
to “lithos” and data paths in general are CRREL/Dartmouth specific. All relevant files are available 
upon request. 
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Figure 5.4.2: Contents of test plot #3 blind test suite #2. 
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October 30, 2005 


Figure 5.4.3: Contents of test plot #4 blind test suite #2. 
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5.4.a Blind Test Suite #2 Results 


The results for the multiobject NSMS inversion are presented on the next two pages. 

The inverted results for this multi-object NSMS method are mixed but encouraging. The it- 
erative multi-object approach here is a slow, somewhat “brute force”ish method. Also remember 
that this data has background subtraction issues as outlined in the last sections. Even with these 
difficulties, the technique here was able to correctly identify individual UXO in about half of the 
cases presented (the green highlighted boxes). For many other cases, the method correctly distin- 
guished between UXO and clutter (yellow boxes), but only rarely produced a full misclassification 
(red boxes). 
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GEM 3D plus processing results for multi objects : 


Blind test # Object #1 / correct object #1 | Object #2 / correct object #2 

1: 

Tp2-1 
2.75 inch / 81mm 
Empty hall 

2 

Tp2-2 Frag CS-10 

3 

Tp2-3 

4 

Tp2-4 


Scrap H-83 / Frag CS-6 


Tp3-1 | | Ϊ 


381 / 105mm HEAT 


Scrape H-259 / 
Frag LG-19 
Or a clutter 


Tp3-2 
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7 
Tp3-3 

Sibert Scrap H-83 / Clutter 

mangle CS-1 
8 
Tp3-4 
9 Signals are too weak to get meaningful information 
Tp3-5 
10 
Tp4-1 
11 
Tp4-2 

S1 sphere 

12 
Tp4-3 

105 HR / 2.75” rocket 
13 ! 
Tp4-4 

BDU-26 / 60mm 
381 /2.75” rocket 
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6 Conclusions and Implications 


SERDP project MM-1537 entitled “Handheld Frequency Domain Vector EMI Sensing for UXO 
Discrimination” is complete. The SERDP project MM1537 has resulted in a unique frequency 
domain instrument that adds transverse receivers to the standard GEM-3 platform while also in- 
corporating several advances in hardware and software. The most notable other hardware advance 
aside from the transverse receiver coils (the 39) part of GEM-39*), a positioning system based 
on the primary field of the instrument itself has been designed and installed which yields subcen- 
timeter accuracy out about 2.5 meters from the beam containing two triaxial receivers. 

Alongside the hardware advances made during this project, the accompanying software and 
models are also of note. The normalized surface magnetic charge (NSMS) model was advanced 
and adapted for the GEM-39*. As well the Standardized Excitations Approach (SEA) was sim- 
ilarly adapted and applied to ]Jgemd data. These physics based rigorous models boast arbitrary 
fidelity and can take into account any primary field and predict the response from any physical tar- 
get. A method to extract an absolute scale which translates GEM-3Y* ppm into physical quantities 
(Amps/meter) was also derived and published. 

The GEM-39%" is now a fully functioning EMI instrument capable of vector sensing of mag- 
netic anomalies while being well located within a limited range. The transverse receivers add 
critical data diversity and information to rigorous models and help in inversion routines. For pre- 
cise, queued interrogation of anomalies, the GEM-39%7 provides diverse, accurate, frequency do- 
main data of the secondary EMI field suitable for inversion and discrimination with high fidelity, 
rigorous models. 
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Α Appendices — Supporting Data 


B List of Scientific/Technical Publications 


Publications produced in whole or part with funds from this project during FY06 are [5, 98, 105, 
110, 111, 121, 126, 134-139]. 


C Other Supporting Materials 


C.1 Data Acquisition Sequence from the GEM-39* 
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Running Gem-3D with GPS unit. 


— Lond, [More than al εἶ 
Ἐπ [meter away ) 


You must follow these instruction to the letter or you will not have good data. 
Turn on: 

1.) Gem-3D -- Battery or Power Supply — 11.6 volts DC 

2.) GPS - Battery or Power Supply — 12 volts DC 


3.) Computer 


4.) Start: WinGem v4 


WinGEMv4 
Press the Transmit button — Then the Null Em Button — zero it. 


“P WInGEM - G3D-1 = ie) x! 
Run Logging Show 


30 90 210 450 1170 3930 13690 3003047010 90030 


eady [ {το eps:eRi ἴση: 123 ἸΕχξιρουν. ἱἘΜ: 000 ἴπ:ῷ (FID: 0 [00:00:11,155 
Bistart| | 4 | 9 Ὁ Εἴ ΤΕ wincem - 650-: | Peasy Button Manager Ϊ MIQXS 9SO EH izszam 


5) Start: GEMPOS ad 
GEMPOS 
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6.) Press the Transmit button — Then the Null Em Button — zero it. 


Pare] 


Run Legging Show. 


Send GPOS data 


54 Jooss21-s213 Alpha: -0.0136231 

»: [azzoe7 Beta: -0.240788 

as festoea PosS: 1,63119e+010 

τ τὶ [Ξτξ ϑββ FQ: 6.40994 wee 
Ready Γ r= fersia [izae [BAT:? ες |Em: coco lim? |FID: 0 [o0:00:20 
Hser|| 4 ΕΞ ΕΞ «Ὁ (2) ΕἾ (Ὁ) | ΕΞ wincem - 630-1 | easy sutton Manager | Eacemeos-tocarone | [PUSIRL@ PEO CFS test am 


dis teas ἐς : : -[πΠ| x 
7.) Then minimize reduce screen by hitting the minus sign. Ἀ ΣΙ 
Do not shut off, to silence beeping. 
8.) Go back to WinGem - you are ready to test - Press the Store Button 
You can use automatic file names or create your own. 
: ParTES 


Run Logging Show 


τῷ Time Graph 


[Al Sensors _~]|| | EFF EEE ONS 0N | 
a Pause ln 
= τα 00:55:45 00:55:50 
Save in: [{{:2 cre! =] « ee 
fe ΕἼ aaa.survey.GEM.IDO1.csv 


7 [ΕἸ bbb.survey.GEM.IDO1.csv 
My Flecent 
Documents 


Ss 


Desktop 


2 


My Documents 
μῳ 


My Computer 


« 
Wabi) File name: 1_gemad.survey.GEM.csv 
Places 
Savesstype: [Survey Files (" Survey.GEM.IDO1.csv] 
23 sor Tro Sas 


ft eps:eri|Sn: 123 |Ext.Pow. ΕΜ: coo jin:? FID: 0 {00:04:15,164 


|\43 ᾿Ξ Ὁ δ J —— 680 1 | Easy Button Manager | Ecempos-tocator2 ἰ πη πὴ μῆς ΤῸ) POD CF 12:56 ΑΜ 


9.) Press Store again when test is finished. 


10.) Then minimize screen by hitting the minus sign. -[Π| x} 
Do not shut off. Ἀ 


11.) | You must run Merge after each test or data will be lost. 


12.) Start: Merge 


merge 
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13.) 


a Untitled - MERGER " 


File Edit View Help 


Press Select Survey — fine the test that you have just finished and open it 


= 16) x! 


Dae] αὶ Ξ5 ΚΞ ἐρ] Φ 


Available Devices/Sensors 
[ Click to discard } 


rE 
Look in: | € 9-23-2008 τὶ 


21x) 
+ © ce EB 


fe 


My Recent 
Documents 


Desktop 
My Documents 
- 
Sr 
<= 
My Computer 


File name: 


δ. 
Peete Files of type: 


Places 


Ready 


14.) 


=) 1_gem3d.survey.GPS.csv 
ΓῚ 2_gem3d.survey.GPS.csv 
=) 3_gem3d.survey.GPS.csv 
(5) 4_gem3d.survey.GPS.csv 
(=) S_gem3d.survey.GPS.csv 
(=) 6_gem3d.survey.GPS.csv 
ΕἾ 7_gem3d.survey.GPS.csv 
=) 8_gem3d.survey.GPS.csv 
ΕἾ 9_gem3d.survey.GPs.csv 
©) 10_gem3d.survey.GPS.csv 
=) 11_gem3d.survey.GPS.csv 
=) 12_gem3d.survey.GPS.csv 
(=) 13_gem3d.survey.GPS.csv 


ΕἸ 14_gem3d.survey.GPS.csv 
ΕἸ 15_gem3d.survey.GPS.csv 
ΕἸ 16_gem3d.survey.GPS.csv 
ΕἸ 17_gem3d.survey.GPS.cs 

ΕἸ 18_gem3d.survey.GPS.csv 


“survey.GPS.csv τὶ 


[Survey Files [*.survey.GPS.csv) 
I Open as read-only 


ΠΞΞ Cf 


Then press Merge LPS-3d Survey button. When Merge is finished it will close 


automatically. 


"Untitled - MERGER 7 - (oI x| 
File Edit View Help 
Π ὦ || ὁ Be|s|? 
Available Devices/Sensors 
{ Click to discard } 
I~ Export to XYZ (Geosoft) [¥ Keep Header Notes Γ Keep Surfer Header 
GPS Settings Smoothing Filter Background Remover —— 
8: fo 

UTM zone: [β΄ ΓΓ Filter before merge ΓΓ Filter before merge τὰ : 

UTM hemis: [N Median Size [7 Median Size [77 Tine Offset ine] Κα 

UTM sys: [wast 984 Mean Size Ρ ἤρα 5:5 Ρ Sensor Offsets From GPS Center Point 

ΠῚ >YLD(-1.05,-0.18) ¥ | 
E:A\DATA\10-41-2008\8_010CT08_124635,SURVEY. ΝΣ 
Survey Type 
@ LPS-3D © Standard Merge LPS-3D Survey | 
Ready =a 2 
᾿Ξ 
oe a 
15.) Check data with Window and = Microsoft 
Office Excell 


Explorer 
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16. ) data should look like this. 


= Ox 


Type a question for help | 


:System=LocalPos 


’ - |B) x 
B J K ° | 
;System=LocalPos Ϊ | | | a 
AS ‘Total SamplesinFile) | | | | | | Ϊ = 
- 21: DataLabel: [00] | | | | Ϊ Ϊ 
| | | 
_ALPHA| BETA ἃ LPSTim LPSmsO  PowerLn!1_30Hz Q1_30Hz 1_90Hz Q1_90Hz [1 2ΊΟΗΖ Q1_210HzI1_450Hz 
0.01 29:59.8 77399787 1.32E+09 3.90E+02 HHH 3.90E+02 sHAEEHE 3.88E+02 eeAHEEEH 3 86E+02_ 


0.01 


30:00.0 77399987 1.31E+09 3.89E+02 #HHHHHT 3.90E+02 HHHHHHH 3.88E+02 HHHHHHHH 3.86E+02_ 
30:00.2 77400187 1.3E+09 3 90E+02 ἘΡΗΒΉΜΗΕ 3. 90E+02 ἘΗΗΜΉΜΗΕ 3. 88:02 HHH 86E+02 


| 7400787. Ἴ3Ε:09 3.89E+02 #eHeRHHH 3.90Ε: | 3.88E+02 AAR 3 86E+02_ 
30:01.0) 77400987, 1.27E+09| 3.89E+02 sHHHHEEH 3.90E+02 seAHHEH 3.88Ε:02 HHRHEHH 3.86E+02 
30:01.2 77401187, 1.25E+09 3.91E+02 #HHHHH 3.90E+02 HHHHHHHH 3.88E+02 sHHHHHH 3. 86Ε 1:02. 
30:01.4 77401387 1.29E+09 3.89E+02 #HHHHHHt 3.90E+02 HHH 3. 88Ε:02 HHH 3.86E+02 


©6909 09/09 G9 a9 cD aD 


Not this. 


ΕἸ Microsoft. Excel SEE 


Type ἃ question forhelp ᾧ 


SEM3D.SURVEY.GEA sv "ἢ fos) 
A B N ο | 
;System=LocalPos = 
45 Total Samples In File 
‘Info Date! 2008 Time: 19: DataLabel: [001 
:Note1 
Ὗ, Z| ALPHA, BETA. ἃ ΕΞ ἰΡΒΤΙπὶ LesmsO_ “PowerLnl1_30Hz Q1_30Hz |1_90Hz Q1_90Hz I1_210Hz Q1_210HzI1_450Hz 


|2.38E+09 3 90Ε:02 #aHHHHHH 3.90E+02 HAE 3 89E+02 βΉΒΒΗΘΘΗ 3.87E+02 
|2.38E+09 3.89E+02 ἀβήθθημβ 3.90E+02 HAHEI 3 89E+02 eeHEEHH 3.87E+02 
2.35E+09 3.91E+02 #HAHHEHE 3.90E+02 HHHHHEHE 3.89E+02 AHHH 3.87E+02 
2.36E+09 3.89E+02 HAHAH 3.90E+02 HAHAH 3 89E+02 HAM 3 87E+02 
(2.34E+09) 3.90E+02| semnraHe 3.90E+02 sete 3. 89Ε:02 aeHHeEHE 3.87E+02 
2.36E+09 3.90E+02 seHHHeHH 3.90E+02 seen Ξ 102 se 3.87E+02 
E 3.87E+02 
1BE+09 3: piste 3.90E+02 seen z 89:02 seteneHe 3 87E+02 
2.35E+09 3.90E+02 #HHHHHH 3 90Ε:02 HHH! 3 89Ε:02 HeHHEHHY 3 .87E+02 
(234E+09 3.90E+02 HAAAKHM 3.90E+02 βάβήβθμᾳ 3 89E+02 HARARE 3.87E+02 
(2.37E+09 3.90E+02 aHHHHHHH 3.90E+02 HAHAH 3 89E+02 βήθθηθθ, 3.87E+02 
|2.37E+09 3.90E+02 saemree 3.90E+02 see 3 .89E+02 seme 3.87E+02 | 


Δ τσ δῶ FANE. AD ἀμ, 2 ANE.AD HuwedauH! FONE Δ’. HuKABBuH 2 ΟἿΣ (ΔΑ 


ΠΕ ΝΣ 


κι κι κι κ el ν κι κν 
δ] «ἱ οἱ οἱ ἡ οἱ οἱ ᾿ἰ οἱ ὦ ὦ οἷς 
ΞΕ Ξ 
νιν. 


17.) You are ready to start your next test. 
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C.2 Beacon Positioning System Construction Details 
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Progress Report (October 2006) 


Progress Report 
on 
Handheld Frequency-Domain Vector EMI Sensor for UXO Discrimination 
Contract Number: W913E5-06-C-0011 


April 2007 


from 
Geophex, Ltd. 
605 Mercury Street 
Raleigh, NC 27603 


1. Progress Summary 


This report covers our progress on the Beacon navigator. CRREL issued a contract 
modification early 2007 so that this task is completed in FY07. As of this reporting, we have 
completed the physical construction of the coils and housing. Early this month, we started 
electronic wiring and assembly. 


2. Navigator Geometry and Construction 


The Beacon Navigator is to be used for accurately locating the GEM-3D active EMI sensor 
Geophex built for CRREL. Figure 1 shows how it works: the GEM transmitter acts as an active 
beacon that can be located by a pair of receiving coils in the vicinity, each of which consists of 
three orthogonal coils. This method accurately produces the sensor's position, as well as the 
sensor's orientation as defined by two tilt angles. Tests have shown that the method produces 
accuracy of better than a centimeter when the receiving coils are at a range of several meters. 


(X1, ΥἹ, 21) 
Two 3-Axis _— 
Receiving Coils “ 
(ΧΟ, YO, ZO) 
(X2, Y2, Z2) α 


gone 
᾿ς Β Transmitter 


Figure 1. A pair of 3D coils measure the magnetic field generated by the EMI sensor. The 
six coil outputs are sufficient to determine both the sensor location and tilt angles. 


We built in 2006 a prototype Beacon Navigator housed in a 10cm x 10cm x 1.5m PVC tube 
that contains a 3D coil at each end. The output of the prototype navigator is routed by hardwire 
(~10m long) to the GEM-3 sensor. The Navigator computes the sensor location using the 
voltage output from the coils and sends the results to the GEM-3 that merges the location data 
with the sensor data. Under this CRREL project, we will refine the existing navigator by 
improving and simplifying its performance, and producing an advanced navigator. We will then 


= SSS 
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integrate the Navigator with the GEM-3D, using the Bluetooth wireless communication so that 
it will be completely detached from the sensor operator. 


Calibration: Since the relative geometry of the two 3D navigator coils is fixed, only a few 
ground-truth points are sufficient for a full calibration. One-time calibration is sufficient for its 
lifetime usage. 


3. Progress Status 
3.1 Navigator Sensing Coils 


Table 1 shows the preconstruction design figures of the three nested coils. The z-axis coil (the 
middle of the nest), for instance, has a mean winding cross-section of 6.8 cm x 6.8 cm and has 
2,000 turns of #34-guage wire. The turns for other nested coils are adjusted to have a similar 
area-turns product to equalize their voltage outputs. Table 1 also shows predicted inductance 
(in Henry) and resistance (in ohm) for each coil. 


As designed 

Rx-Z Middle Section 
2000 turns 0.873 AW G#34 9.25 area x turn (m*m) 
0.068 a (m) 0.334 inductance (H) 64 x-sec (mm**2) 
0.068 b (m) 474.9 resistance (ohm) 0.703 time constant (msec) 


Rx-X  Outermost 


1784 turns 0.873 AW G#34 9.25 area x turn (m*m) 
0.072 a(m) 0.290 inductance (H) 57.088 x-sec (mm**2) 
0.072 b(m) 448.5 resistance (ohm) 0.647 time constant (msec) 


Rx-Y — Innermost 


2260 turns 0.873 AW G#34 9.26 area x turn (m*m) 
0.064 a(m) 0.387 inductance (H) 0.766 time constant (msec) 
0.064 b (m) 505.1 resistance (ohm) 


Table 1. Design figures of the three nested coils 


Figure 2 shows actual coils constructed based on the design figures. Each of the two sensing 
coils is wound on a light-density, 3-inch plastic cube that contains a pair of recessed grooves 
for coil on each axis. The coils are to be mounted at the ends of a 3” x 3” x 5 ft square PVC 
tube, as shown in Figure 3. 


Table 2 shows the measured inductance, resistance, and resonance frequency of actual coils 


that were built based on the design figures. The predicted and actual are close within 
manufacturing accuracy. 
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4 


Figure 2. Two 3-D nested sensing coils with preamps attached 


Figure 3. Square PVC tube housing 


es 
half half total half half total 
0.116 0.117 0.323] 0.117 0.117 0.334 ]ind (H) 
240 241 481 235 235 470]}res (Ohm) 


res fea (2) 


a ees ee 


0.102 0.101 0.288] 0.101 0.101 0.287 ind (H) 
a ee ΕΝ 
137 137 0.393] 0.138 0.138 0.394 ]ind (H) 


Table 2. Measured inductance, resistance, and resonance frequency of the actual coils 
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3.2 Sensing Electronics 
The final electronic design for the navigator is still in progress. Two primary design issues are: 
e Wireless communication with the GEM-3D and 
e Automatic gain control scheme for each preamp to prevent input saturation when the 
GEM- 3 is close to the navigator. 
The second issue is a newly-discovered: when the GEM-3D is at less than, say, 1 m from the 
navigator, the front-end voltage becomes high to saturate the preamp. We are currently 


working on an automatic gain-control scheme to deal with this issue. 


Figure 4 shows the main DSP section that will process the 6-channel time-series data. 


Figure 4. DSP electronics for the navigator 


4. Plan for the Next Reporting Period 


We will continue the electronic assembly of the navigator. We expect to fire it up within the 
next month and to begin the test and calibration phase. 


5. Estimated Percentage of Contract Completion for Year 2 Task: 40% 
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C.3. GEM-39* Data Processing - Static and Dynamic 
C.3.a_ Introduction 


This short document explains the internal structure of the program for data read-in and analysis of 
GEM-39" files. 
Two cases are considered: 


1. The “Beacon sliding on a grid” data taken in the lab by Ben Barrowes on 10/7/08. This set is 
stored on lithos in directory 


DATA/gem3dplus/3D_10-7-2008 _Ben/ 


and is henceforth called Ben. 


2. The “Beacon free form” data taken over the test plots by Kevin O’Neill on 10/23/08. This 
set is stored in 


DATA/GEM-3D+_OCT-23-08_TESTPLOTS EXPORT/ 


and is called Kevin from now on. 


From now on when I refer to a “file” I really mean a collection of three .csv files: 1001, 
ID02, and 1003, corresponding to the GEM-3Y*’s three receivers. 


C.3.b Reading the data 


1. All filenames start with an identifying number. The ranges (of usable files) are [21 : 49] 
(29 files) for Ben and [1,3,4,8,9,12 : 18,21,22] (14 files) for Kevin. The variable Ncase 
ranges from 1 to 29 (14). 


2. Some files have a strange feature: ID01 has one more line than either ID02 or 1003. This 
happens both in Ben’s and in Kevin’s data. 


Kevin’s writeup mentions a “slight timing asynchronization” in the data he studied ( (2-1) , 
or 3. κ in our present notation), which happens to be one of those cases. He suggests two 
solutions: interpolating or just ignoring the small differences. 


I have taken the second route and simply deleted the first row when the set is problematic. 
At some point I tried interpolating (since each datafile includes the time at which each point 
was measured) but did not see much of a difference. I can redo this if necessary. 


3. The previous version of the data-reading procedure made some use of UNIX system calls 
to determine the filenames and read in the values. I fixed that feature, which made it both 
system-dependent and slow, and replaced it with Matlab-only functions. 


As it stands, the procedure can read . csv-like files that need not be all numeric. The only ad 
hocs here respond to Kevin’s having deleted most of the headers and the “LPSTime” column 
that presumably was preventing him from using csvread. 


4. All three files corresponding to a given “file” are read at once and stored in a 3D array. 
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5. Each receiver provides its own position/tilt information. I have kept this as both a 2D matrix 
(e.g., Xpos, with leading capital, is Not x 3, where Npt is the number of measurements) 
and as an averaged array (xpos, all lowercase). 


We do note that question arises as to the accuracy of the measurement locations given by the 
beacon system. A look at Xpos and the other two such variables reveals two columns that are 
identical and distinct from the third. The differences are most notable for the x-coordinate, 
where they reach 3 cm and are not random but periodic and roughly 2/2 out of phase with 
the actual x-values. (See Fig. C.3.2 for a plot.) 


6. We also delete rogue points in which the location information was patently wrong, for ex- 
ample featuring large discontinuities. In Ben’s case, it is straightforward to find those points, 
because they had a yaw angle of exactly zero. In Kevin’s case we had to find them by hand. 


C.3.c Background subtraction 


C.3.c.(1) Initial approach We start by noting that we neglected the angular information in the 
background subtraction. In all cases some effort was put into holding the sensor as level and 
well-aligned with the beacon coordinate system as possible. 

Ben took a whole background file, Ncase = 2, which makes the subtraction straightforward: 
just read in that file, average over all points, and subtract the average from each data file. We do 
not have to be too observant of the height in this case, as there is no soil. 

Kevin’s files are a little different, since he takes the background values from the data them- 
selves. (Also, his “calibration” file is over a shotput.) In his words, 


The program next asks you to enter some data point number around which it will aver- 
age (+2 points) to obtain a background to subtract from all data treated subsequently. 
It’s often hard to come up with a good choice, which has implications for deeper pro- 
cessing. Here we’re just doing some exploratory visualization so it’s not so important. 


Now, however, is the time to perform the “deeper processing” he mentions. We cannot use Ben’s 
background file, since that does not include soil response. 
On the other hand, Kevin chooses his background point thus: 


In any case, in the favorable case shown here, point #25 occurs after the surveying has 
gotten under way and is around the beginning of the first sweep; and there isn’t much 
going on in either of the two response measures. Thus we enter this and the program 
averages over points 23—27, plotting the resulting background profile so you can see if 
anything wacky is going on, i.e. more strange than usual. 


This gives rise to two related possibilities for finding the background: 


1. Decide upon a fixed number of background points, sort the measurements in decreasing order 
of distance from the center, and use the (averaged) furthest measurements as background. 


2. Carry out a multi-file survey and use the very furthest measurements as background, after 
averaging them. 


Now, how do we determine an “optimal” number of background points that will not skew the 
measurements? We must look at the way the (frequency by frequency) average for the furthest 
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points change as we increase or decrease Nbk; also, we want to see how it compares to the field 
values themselves. 

It turns out the background values so chosen exceed the combined background + UXO mea- 
surements in a few cases, and moreover the response is not flat. This may well be due to the height 
dependence of the soil response, which is nominally ἢ. On the other hand, the response of a 
dipole is approximately d~°, which lends confidence to our taking the points presumably furthest 
from the target. 

The user can control (dictate) Nok, the number of points to take as background. The program 
finds a “center” point defined by 7, = χοῦ" y-¥ + 202, where x, and y, are the averages of x and y 
and zg is the lowest possible value of z, centers the positions there, and measures the radial distance 
from that origin to every measurement point. It then counts the Nbk points furthest from the origin, 
averages their field values, and uses this as a background. The function takes Nok = 20 as default, 
though that value can be changed. 


C.3.c.(2) Two-height approach The method described above has some flaws. First, it mashes 
together “background” data corresponding to different heights, which can bias the results. Also, 
the distance from the “‘center,’ which we estimate as the average of the x- and y-locations, may 
not—and in multi-object cases will not—coincide with the center of anything, and thus it may turn 
out that the background points are in fact right on top of a target. 

We then had to find a background-subtraction technique that would try to overcome these prob- 
lems. In the end we implemented an ad hoc two-height approach. 

Kevin’s measurements consist of two zigzags, one close to the ground making five passes over 
the measurement area and the other retracing the original in the opposite direction and some 10 cm 
above. The measurements were taken by hand, so the separations between the two zigzags are 
not necessarily clear-cut, and in some cases there are “higher” points that are actually closer to 
the ground than some in the “lower” zigzag. In our procedure, however, we assume two distinct, 
consecutive heights. 

We start by finding by hand a single point, halfway in the transition region between the two 
zigzags, and divide each data set into the “before” and “after” intervals. We then find a five- 
consecutive-point interval for each of the zigzags that represents the best compromise between the 
following desirable (and often conflicting) features: 


1. The interval is far from the center of the measurement trajectory. 
2. The points lie at roughly the same height, one representative of the corresponding zigzag 


3. The five points lie along a fairly smooth path so the mean of the field values is a reasonable 
average. 


4. The resulting fields will have these characteristics: 


e The z-field should be positive in as many points as possible. 


e The x- and y-fields should have roughly as many positive values as negative, and the 
zero crossings should be roughly at the center of each sweep. 


We concentrated on the quadrature parts of the fields, since they are presumably less affected 
by ground effects and did not have the large anomalous dc shift that the inphase parts had. 
Moreover, the quadrature response is positive at all frequencies, so any negative amplitudes 
will be due to geometry. 
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5. The above rules should hold in as large a range of frequencies as possible. 


It is this method that we used in the end to prepare the data for further processing. 
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Normalized center position 


y -0.5 -0.6 x 


Figure C.3.1: Sensor position in the example data run. The convention in this and the following figures is: 
the sensor “warms up” as it takes data: i.e., the trajectory is initially blue and transitions into red during the 
run. The background points are black dots. 


C.3.d An example of the first approach 


As an example we include our analysis of Kevin’s file Ncase = 2, corresponding to location 2-1. 

For this example we use the first approach to background subtraction. We choose Nok = 10 
in this case; the points happen to lie at the very end of the data run, and appear as black dots on 
Fig. C.3.1. The black star in that plot is the center from which the background points are most 
distant. 

Figures C.3.3 and C.3.4 respectively show the inphase and quadrature parts of the background- 
subtracted secondary field at f = 450 Hz. On the left-hand-side panel we have the raw values. The 
dashed lines represent the background value, computed by averaging the field values component 
by component at the Nbk points, which appear in the figure as dots. On the right-hand-side panels 
we have the values obtained after subtracting the background. 

The quadrature parts are further depicted as “profiles” on Fig. C.3.5. 
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Figure C.3.2: Discrepancy in the x position between the first and third data files of a set. 
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Figure C.3.3: Background-subtracted field at f = 450 Hz, inphase part 
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Figure C.3.4: Background-subtracted field at f = 450 Hz, quadrature part 
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Figure C.3.5: Background-subtracted 3D fields at f = 450 Hz. 
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Figure C.3.6: Quiver plot showing the projection on the x-y plane of the quadrature part of the secondary 
field at f = 450 Hz. Joining the tails of the arrows can give a approximate estimate of the locations of the 
targets. 
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Normalized center position 


Figure C.3.7: Sensor position in the second example data run. 


C.3.e Example of the second approach 


To illustrate the second, final background-subtraction method we take another example, location 
2-4, 

Figure. C.3.7 again plots the locations, with labels and hollow markers identifying every 10th 
point. The background points can be seen as solid black dots, to the left and right of the figure. 

Figures C.3.8 and C.3.9 respectively show the inphase and quadrature parts of the background- 
subtracted secondary field at 3930 Hz. The description is as before: the left-hand-side panel shows 
raw values, the dashed lines represent the (two-level) background level, the dots depict the points 
used to find the background. On the right-hand-side panels we have “normalized” fields. 

The field profiles of Fig. C.3.10 look two-peaked, and the impression that this plot may have 
two objects is strengthened by looking at the quiver plot in Fig. C.3.11. 
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Figure C.3.8: Background-subtracted field at f = 3930 Hz, inphase part. 
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Figure C.3.9: Background-subtracted field at f = 3930 Hz, quadrature part. 
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Figure C.3.10: Background-subtracted 3D fields at f = 3930 Hz. 
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Figure C.3.11: Quiver plot showing the magnitudes of the x- and y-components of the quadrature part of 
the secondary field at f = 3930 Hz. Joining the tails of the arrows can give a approximate estimate of the 
locations of the targets. 
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